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ABSTRACT 


A relatively simple time domain method is developed to calculate the time of arrival for 
radar signals The error present in the estimate of the time of arrival for a single pulse 
and a burst of pulses are developed and the effects of SNR, PRF. pulsewidth. and 
sampling frequency are examined. T -ival is used with multiple sensors and the 

Kalman filter to estimate the location of the em'\er. Algorithms estimate the location of 
an emitter given the Time Difference of Arnval ( r^X)A) of a single pulse as well as the 
TDOA of bursts of pulses received as the emitter scans past the a-ceivers. The algorithms 
were tested on simulated data. 
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I. INTRODUCTION 


A. BACKGROUND 

During the Persian Gulf War UAVs (Unmanned Aerial Vehicles) proved useful by 
providing battlefield surveillance, targeting information, and naval gunfire spotting. 

They were able to provide the battlefield commander with up to the minute intelligence 
and offered the battlefield commander an over the horizon intelligence gathering asset 
that possessed the flexibility that he required to respond to the constantly changing 
circumstances and priorities of the battlefield. UAVs also provided a medium range 
aerial reconnaissance and intelligence gathering capability that was not hindered by the 
long lead time required for other reconnaissance assets. Current UAV capabilities are 
primarily limited to optical and infra-red surveillance, but the UAV would be very useful 
in electronic intelligence (ELINT) gathering. UAVs could remain on station listening to 
the enemy's electronic emissions for long periods of time, deep in enemy territory, 
without risking aircraft and pilots. 

One application of electronic intelligence gathering is locating radar emitters from 
time difference of arrival (TDOA) observations. A time domain technique is presented 
here that uses the Kalman filter to estimate the location of a radar emitter from the time 
difference of arrival for a burst of pulses between two or more receivers, and time 
difference of arrival of individual pulses between two or more sensors. The background 
information and models used to develop the Kalman filter are presented, and the resulting 
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Kalman filter algorithms are tested for a variety of problem configurations. Conclusions 
are drawn about the usefulness of the algorithms and recommendations are made for 
further testing and evaluation. 

B. RADAR EMITTER SIGNAL CHARACTERISTICS 

The signal characteristics most important in the emitter location problem are the 
Pulse Repetition Interval (PRI), and the pattern and rate at which the emitter scans a 
target. These characteristics determine which estimation technique will be most 
effective, and they affect the ability of these algorithms to accurately estimate the location 
of the emitter. 

The PRI is the period of time between successive pulses. This value can range from 
tens of milliseconds for long range search radars to microseconds for pulse Doppler 
radars. Range in many radar systems is determined by measuring the time required for a 
pulse to strike a target and return, and the unambiguous range of a radar system is 
determined by the length of the PRI. The longer the PRI the longer the unambiguous 
range of the radar. Similarly, the PRI affects the allowable separation of the sensors used 
to detect the pulse TDOA. The longer the PRI the greater the allowable separation of the 
sensors. 

The vast majority of radar systems are designed to scan a volume of space to look for 
targets. All of these systems use a search mode that is used to scan the volume of space 
around the emitter in a regular pattern. Most radar emitters spend the majority of their 
operation in this mode and it is in this mode that most emitters can be passively detected. 
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The method that the emitter uses to scan and search for a target is imp>ortant in the 
problem of estimating its location from its emitted signals. Radar systems typically scan 
their main beam mechanically or electronically. 

Mechanically scarmed radar systems designed for target location typically scan in 
azimuth in a circular pattern. Their antennas are rotated by motors and actuators in a 
circular manner, and are designed to provide 360 degree search coverage. Special radar 
systems designed for altitude finding, or tracking and fire control will not necessarily 
scan in azimuth alone. Most of these systems use electronically scanned antennas to steer 
the main beam of the radar, and can scan in elevation as well as azimuth. These systems 
can have very complex scan patterns that can be changed as required. This thesis will be 
primarily concerned with mechanically scanned radar emitters that scan circularly in 
azimuth. The majority of the emitters encountered will be of this type. The principles 
developed to locate these emitters can also be applied to emitters with more complex scan 
patterns. 

As an emitter scans its main beam past a target, the amplitude of the pulses that 
strike the target will generally vary from zero to a maximum that occurs near the center of 
the main beam. Figure 1 shows the typical burst of pulses detected as the emitter scans 
past a receiver. The frequency of these bursts is a function of the scan rate of the emitter, 
and the width of the burst is a ftmction of the scan rate of the emitter and the beamwidth. 
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Figure 1. Burst of pulses from a scanning emitter 

C. EMITTER LOCATION PROBLEM DESCRIPTIONS 

Two basic approaches are used to estimate the location of the emitter from the 
TDOA of its signal. In the first approach the TDOA for the burst of pulses detected as 
the emitter scans past receivers is used to estimate the location of the emitter. The second 
approach uses the TDOA for individual pulses between two sensors on widely separated 
platforms to estimate the location of the emitter. 

1. Burst TDOA Problem Description 

The first TDOA filtering problem assumes that multiple receivers are used and 
that they are widely separated. The TDOA for the burst of pulses detected at each 
receiver are measiued and used to estimate the location of the emitter. In this burst 

TDOA filtering problem, the following assumptions are made; 

1. Enough information is obtained fix>m the received signals to associate TDOA 
observations with the proper emitter, and thus the multi-emitter problem is 
reduced to a single emitter problem. 







2 The receivers are referenced to a common ume base All time of amvai 
measurements are referenced to this common time base 

3 GlohaJ Positioning System (GPS) is used to determine the locations of the 
receivers. These estimates of the receiver locations are assumed to be exact and 
errm introduced by the GPS is ignored. 

4. The emitter scans in azimuth at a constant rate. 

5. The PRl is constant for the entire length of the burst. 

The receivers may be mounted on UAVs, aircraft, ground sites, or space platforms. For 
the UAV mounted receivers, the desired maximum detection range is assumed to be 
500 kilometers. The algorithms developed are tested for ranges of 500 kilometers as well 
as shorter ranges of 150 kilometers and 30 kilometers. 

2. Single Pulse TDOA Problem Description 

The pulse TDOA problem estimates the location of the emitter by measuring the 
TDOA for an individual pulse between two sensors. The sensors are spaced close enough 
so that the TDOA can be measured for a single pulse without ambiguity. The following 

assumptions are made in the problem development: 

1. Enough information can be obtained from the received signals to associate TDOA 
observations with the proper emitter, and thus the multi-emitter problem is 
reduced to a single emitter problem. 

2. The receiver platforms are referenced to a common time base. All time of arrival 
measurements are referenced to this common time base. 

3. Each receiver platform has two sensors and the TDOA for individual pulses 
between the two sensors can be measured accurately. 

4. GPS is used to determine the locations of the receiver platforms. These estimates 
of the reiver locations are assumed to be exact and error introduced by the GPS 
is ignored. 

The algorithms are tested for ranges of 500 kilometers as well as shorter ranges of 150 
kilometers and 30 kilometers. 
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II. TDOA ESTIMATION 


A. TIME OF ARRIVAL ESTIMATION AND ERROR MODELS 

The Kalman filter can be used to estimate the location of the emitter from TDOA 
observations. Kalman filter theory, however, assumes that the noise corrupting 
observations is Gaussian distributed. Before applying the Kalman filter to the emi tier 
location problem, an analysis of the error statistics of the TDOA observations is 
necessary. This determines if the noise present in TDOA observations can be 
approximated as Gaussian. One method of estimating the TDOA for emitter bursts and 
pulses is presented. The effect of noise on this estimation method's accuracy is explored, 
and the probability densities calculated are compared to Gaussian distributions. 

1. Time of Arrival Measurement 

TDOA measurements between widely separated receivers and sensors are used to 
estimate the location of the emitter. Since the receivers need to be widely separated to 
improve the accuracy of the location estimates and the orthogonality of the TDOA 
observations, physically linking the receivers together is not feasible. Therefore, it is 
important that the receivers are accurately referenced to some external time base. An 
onboard time base can be used to estimate the time of arrival (TOA) of signals, but drift 
would introduce error into the estimates of the TOA and subsequently the TDOA 
observations. To minimize drift, internal clocks aboard the receivers can be periodically 
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corrected with references to a more accurate external time base. If the receivers are 
mounted on UAVs equipped with GPS, the GPS system time could be used as the 
external time base. This approach is assumed for the balance of the development. 

The TDOA observations in the burst TDOA filtering problem are on the order of 
milliseconds and tens of milliseconds. The error typically present in the GPS system 
time is much smaller than the typical TDOA observations and synchronizing all of the 
receivers to the GPS system time would not introduce appreciable error into the 
observations. For these observations, the receivers are assumed individually 
synchronized to the GPS system clock. 

For the pulse TDOA observations, the length of the typical observations are a 
few microseconds or less. Using the GPS system time to estimate the time of arrival of a 
pulse at a sensor would introduce unacceptable error in the TDOA observations. This is 
because the error present in the GPS system clock would be on the same order as the 
observation it is being used to measure. To avoid ambiguity in the pulse TDOA 
observations, the sensors must be placed relatively close together. For the balance of the 
developments here the sensors will be assumed physically linked to the receiver platform, 
and the TOA observations for both sensors will be referenced to the same on board time 
base. 

2. Pulse TOA Estimation 

Envelope detection of the pulses received by the sensors is assumed. The 
envelopes of the received pulses are sampled and digitally encoded. An estimate of the 
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TOA for a pulse can be obtained from the time at which a sample is detected above a set 
threshold. This method of estimating the TOA is highly sensitive to the sampling rate 
and noise present in the sampled pulse. Spurious noise that is above threshold gives false 
TOA indications. Algorithms can be developed to test for false indications by examining 
multiple samples and announcing the arrival of a pulse only if a specified number of 
above threshold samples are detected. A more accurate estimate of the TOA can be 
obtained by integrating the samples of the pulse and calculating the TOA of the centroid 
of the pulse. The improvements gained from integration are similar to the improvements 
obtained in radar systems when pulse integration is used. All of the information present 
in the pulse is utilized and the effective signal to noise ratio is increased. 

The time of arrival of the centroid of a sampled pulse can be found from 
equation (1): 


TOA = 


£tOc)z(k) 

ks\ 



( 1 ) 


Where TOA = the Time Of Arrival (TOA) of the pulse centroid, 

z(k) = the magnitude of the kth sample, 
t(k) = the time at which the kth sample was taken, and 
N = the number of samples in the pulse. 

In the following development of the error present in the estimate of the pulse 


TOA, detected pulses are assumed to have the shape shown Figure 2. To these pulses. 


white, Gaussian distributed noise is added and its effect on the TOA is calculated. 
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Figure 2. Pulse envelope and sampled pulse 


The sampling interval is assumed to be much smaller than the pulse width so a 
sufficient number of samples of the pulse are taken and the sampled pulse is a reasonable 
representation of the original pulse. The sampling is assumed to be triggered by, and 
aligned with, the rising edge of the pulse. This assumption ignores the error introduced in 
the TOA when the first samples of the pulse are taken at some random point within the 
pulse. 
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Each sample of the pulse is considered a random variable as shown in 
equation (2). 

2(k) = Zo(k) + v(k) (2) 

Where z(k) = the random variable representing the amplitude of the kth sample, 

Zo(k) = the deterministic variable that is the true amplitude of the kth 
sample, and 

v(k) = the white Gaussian random variable that represents the noise 
present in the amplitude of the kth sample. The noise has zero 
mean and a variance of V, i.e.. E[v(k)^] * V. 

Since each sample of the pulse is modeled as a random variable, the time of arrival of the 
centroid of the pulse, calculated in equation (1), is also a random variable. A detailed 
derivation of the probability density of the pulse TOA is presented in Appendix A. 

The peak signal to noise ratio is used as a means of representing the noise power 
present in the samples of the pulse. As shown in Figure 2, the peak amplitude of the 
pulse is assumed to be one. The peak signal to noise power ratio is calculated from the 
following equations: 

f Speik "! _ E[Zo(k)^] 

^ N J E[v(k)2] ^ 

Where z^Ck) - the true amplitude of the peak sample of the pulse, and 

v(k) = the zero mean white noise sample with variance V, and 
E[ ] = the expectation operator. 

In decibels the peak signal to noise power ratio is given as: 
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Since the peak amplitude of the pulse is one and the noise v(k) is a random variable with 
zero mean and a variance V, the peak signal to noise power ratio reduces to: 

, 5 , 

The MATLAB program Puldist.m, listed in Appendix E, calculates the 
probability densitiy functions plotted in Figure 3. These curves depict the probability that 
the TOA of a pulse, calculated using equation (1), will fall within a specified range. The 
mean and standard deviation of the probability density fimctions are also calculated and 
Gaussian distributed density fimctions having the same mean and standard deviation are 
superimposed on the plots. 



Figure 3. Pulse TOA probability densities vs. Gaussian 
approximations, pulsewidth 1.0 psec 
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The Gaussian density functions are nearly indistinguishable from the actual 
probability density functions, so the pulse TOA can be considered Gaussian distributed in 
the Kalman filter developments. In Appendix A, a complete analysis is conducted of the 
effects of peak signal to noise ratio and sampling rate on the probability density function 
of the pulse TOA. 

3. Burst TOA Estimation 

As the main beam of an emitter sweeps past a receiver, the receiver detects a 
burst of pulses. The number of pulses received is a function of the beam width and scan 
rate of the emitter. The amplitude of the pulses change as the main beam of the emitter 
sweeps past the receiver. A burst typical for circularly scanning emitters is shown in 
Figure 1. The envelope of the burst can be approximated with the upper half of a sine 
wave. 

x(t) = sm[^®t] (6) 

Where x(t) = the amplitude of the envelope of the burst, 

BW = the beamwidth of the main beam of the emitter in radians, 

SR = the scan rate of the emitter in radians/sec, and 
t = the time in seconds. 

An estimate of the burst TOA can be obtained from the time of arrival for the 
first pulse of a burst, but this method of estimating the TOA assumes that for each 
receiver the same pulse in the burst will be detected first. Because of lower signal 
strength, more distant receivers may not detect some pulses as the start of the burst and 
the estimate of the TOA will be in error by some multiple of the PRI. Estimating the 
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centroid of the envelope of the burst and its time of arrival yields a better estimate of the 
TOA because the shape of the envelope detected by all of the receivers is the same, and 
only the amplitude of the envelope changes as a function of distance from the emitter. 

An estimate of the TOA of the burst envelope is found by integrating the sampled 
pulses and their time of arrival over the length of the burst. If enough pulses are present 
in the burst, the centroid of the pulses lies sufficiently close to the centroid of the 
envelope. Integration of the pulses also reduces the effect of the noise by averaging its 
effect over a longer period of time, and effectively improving the signal to noise ratio by 
ignoring the interpulse periods where only noise is present. 

Since the burst is considered a collection of sampled pulses, equation (1) is used 
to estimate the time of arrival of the centroid of the burst. The equations and algorithms 
used to calculate the probability density function, mean, and standard deviation of the 
pulse TOA are used to calculate these properties for the burst TOA. 

The MATLAB program Burdist.m, listed in appendix E, generates the burst in 
Figiu^ 4 and calculates its probability density for peak signal to noise ratios of IS and 
25 dB. 

The burst has the following characteristics: 

Burst Length: 2.0 milliseconds 

Pulsewidth 1.00 microseconds. 

Pulse Repetition Frequency 6000 Hz, 

Maximum Amplitude 1.00 


13 



Figure 4. A burst of pulses with a 6.0 kHz PRF and 
a 1.0 microsecond pulsewidth 


The envelope of the burst is calculated using equation (6), and the peak signal to 
noise ratios are calculated using equations (4) and (5). The burst in Figure 4 assumes that 
the pulses are spaced evenly throughout the envelope of the burst. If this is the case, the 
centroid of the pulses coincide with the true centroid of the envelope. For emitters with 
long PRIs where few pulses are present in the burst, this assumption may not be valid, 
and the error introduced could be significant. In all further developments emitter 
characteristics are chosen so that this error is insignificant. 

For the probability density function of the burst shown in Figure 4, the mean and 
standard deviation are calculated, and Gaussian distributions having the same mean and 
standard deviation are simultaneously plotted. Figure 5 is a plot of the burst TOA 
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probability density function and the Gaussian approximation for a peak signal to noise 
ratio of 15 dB. Figure 6 is the plot of the burst TOA probability density function and 
Gaussian approximation for a peak signal to noise ratio of 25 dB. 



Figure 5, Burst TOA probability density vs. Gaussian approximation, 
2.0 millisecond burst, 6.0 kHz PRF and 1.0 microsecond pulsewidth 
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Figure 6. Burst TOA probability density’ vs. Gaussian approximation, 

2.0 millisecond burst, 6,0 kHz PRF, and 1.0 microsecond pulsewidth 

The burst TOA probability density functions are nearly indistinguishable from 
the Gaussian probability density' function. The burst TOA can be considered Gaussian 
distributed in the Kalman filter development. In Appendix A a more complete analysis is 
conducted of the effects of peak signal to noise ratio, pulse width, and PRF on the 
probability density functions of the burst TOA. 


B. BURST AND PULSE TDOA 

As demonstrated in the previous sections, the TOA of pulse or burst can be 
approximated as a Gaussian random variable. If the TOA for each burst or pulse is 
considered independent, the TDOA can be calculated using the following equation; 
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TDOA = Za - Zb 


(7) 


Where 2^= TOA of burst or pulse A, and 
Zb = TOA ofburst or pulse B. 

Time of arrival A and B can be bursts generated by an emitter scanning two different 
receivers, or the bursts generated by different scans of the emitter past the same receiver. 
Time of arrival A and B can also be the time of arrival for two different pulses at the 
same sensor or the time of arrival of a single pulse at two different sensors. 

If the TOA is expressed as the sum of a constant value uliich corresponds to the true 
TOA, and a zero mean Gaussian random variable, then the TOA is: 


Za=ZoA+Va 


( 8 ) 


Where Za= the observed time of arrival, 

ZoA= the true time of arrival, and 

Va = the noise present in the observed time of arrival. 

This noise has zero mean and variance Va. 

The noise term Va is considered to be white noise, with zero mean and uncorrelated 
increments. Therefore, the expected value, E[VaVb], is zero. The mean of the estimate of 
the TDOA are: 


Ptdoa = E[za-zb ] = E[Zoa-Zob ] + E[va ]-E[vb ] (9) 


PTDOA = ZoA - ZoB 


( 10 ) 


Therefore, the expected value of the TDOA is the difference between the true times of 
arrival for the two bursts, i.e. the true TDOA. 

The variance of the estimate of the TDOA is: 
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^woA ~ E[ (TDOA - ^tdoa)^ ] 


(11) 


<y™A=E[(VA-VB)M (12) 

Since the increments of the noise are independent, and E[vaVb] = 0, the resulting variance 
of the TDOA estimate is: 

o™A = E[vi + v|] = V* + VB (13) 

Therefore, the TDOA of a burst or a pulse is a Gaussian random variable with a 
mean value equal to the true TDOA, and a variance equal to the sum of the variances of 
the burst or pulse time of arrival estimates. 
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ni. EMITTER LOCATION 


The TDOA for a radar signal is the difference in the amount of time that a radar 
signal takes to reach two receivers. This research utilizes the time difference of arrival 
for bursts of pulses between two widely separated receivers and the TDOA for individual 
pulses between two closely spaced sensors to estimate the location of an emitter. Both of 
these TDOAs are iunctions of the emitter and receiver locations and the emitter 
characteristics. In the approaches pursued here only the two dimensional problem is 
considered, and the effect of altitude is ignored. 

A. BURST TIME DIFFERENCE OF ARRIVAL (TDOA) 

1. Problem Geometry 

The TDOA of a radar signal burst for two widely separated receivers is primarily 
due to the amount of time required for the emitter to scan the two receivers. If the 
assumption is made that the emitter scans only in azimuth, then the TDOA is only a 
function of the scan rate of the emitter and the angle formed between the emitter and the 
two receivers. The scan rate is the angular rate of rotation of the emitter, and the angle 
formed by the emitter and the two receivers can be determined from the coordinates of 
the emitter and receivers. If the scan rate of the emitter is assiuned constant and known 
a priori or from some other estimation method, then the TDOA becomes a function of 
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only the location of the emitter and the receivers. Figure 7 is the geometric relationship 
between the receivers, and the emitter. 



Figure?. Burst TOO A geometry 


Since the scan rate is known, the TOOA for an emitter scanning two receivers is 
only function of the angle 6 formed by the emitters and the receivers. This angle can be 
found by applying vector geometry. If the positions of the receivers are known, the cross 
product of the position vectors from the two receivers to the target can be used to express 
the angle formed by the two receivers and the emitter in terms of the unknown quantities, 
the X and y coordinates of the emitter. The angle 0 is calculated from the following 
equations: 
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RaXRb = iRallRbl sin(0) 


(14) 


M A A 

Ra = (xt - xa)i + (yt - ya)j 


(15) 


Rb = (xt - xb)i + (yt - yb)j 


(16) 


RaXRb = 


A A A 

i j k 
(xt-xa) (yt-ya) 0 
(xt-xb) (yt-yb) 0 


(17) 


sin(0) = (xt-xa)(yt-yb)-(yt-yaXxt-xb) 

[(xt - xa)^ + (yt - ya)^l ‘^[(xt - xb)^ + (yt - yb)^] 

Where Ra = The position vector from receiver A to the emitter, 
Rb = the position vector from receiver B to the emitter, 

0 = the angle formed by the receivers and the emitter, 
xt, yt = the X and y coordinate of the emitter, 

xa, ya = the x and y coordinate of receiver A, 

xb, yb * the X and y coordinate of receiver B, and 
SR == the scan rate of the antenna in rad/sec. 

If 0 is small, the TDOA is found from the following equation. 
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TDOA = 

SR SR 


(xt-xa)(yt-yb)-(yt-yaXxt-xb) 

L [(xt - xa)^ + (yt - ya)^J ''^[(xt - xb)^ + (yt - yb)^] J 


(19) 


2. Loci of Constant TDOA 

From 'quation (19) the burst TDOA for two widely separated receivers being 
scanned by a constantly rotating emitter is directly proportional to the angle formed by 
the two receivers and the emitter. In Appendix B an approximate relationship is 
developed to plot the loci of points where an emitter could lie and produce a specific 
TDOA observation. These loci of constant TDOA are shown in Figure 8. 
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Figure 8. Loci of constant TDOA for an emitter scanning at a constant rate, 
scan rate 360 degrees/sec, receiver separation 1.0 kilometers 

These loci demonstrate that for a single TDOA observation there exist an infinite 
munber of possible location for the emitter. Another observation or measurement is 
required to uniquely estimate the location of the emitter. 

If more than two receivers are used, a TDOA observation is available for each 
combination of receiver pairs. Plotting the loci of constant TDOA for each of the 
receiver pairs, multiple intersections occur, but all of the loci intersect at a common point, 
the location of the emitter. This is shown in Figure 9. 
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Figure 9. Loci of constant TDOA for an emitter scanning at a constant rate, 
scan rate 360 degrees/second, receiver locations as shown 

3. Orthogonality of TDOA Observations 

The performance of the Kalman filter varies based upon the quality of the 
observations available. Generally, the observations with the highest orthogonality 
perform the best. If data processing time or capability is limited, it is advantageous to use 
the observations with the highest orthogonality. The loci of constant TDOA are one way 
of visualizing and calculating the orthogonality of the observations. 

A measure of the orthogonality of the TDOA observations is obtained from the 
dot product of the unit vectors tangent to the loci of constant TDOA. These vectors are 
calculated from the slope of the line tangent to the loci at the estimated location of the 
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emitter. Once the unit vectors tangent to the loci are found, their dot product yields the 
cosine of the angle formed by the two loci. The equation for the dot product is given 
below: 


Rti»Rt2=cos0 (20) 

Where Rt, = The unit vector tangent to the loci munber 1, 

Rt} = the unit vector tangent to loci number 2, and 
0 = the angle formed by the two loci. 

The cosine of the angle formed by the loci is an indication of the orthogonality of 
the loci. If the cosine of 0 is close to one then the loci are nearly parallel, and if the 
cosine of 0 is zero the loci are orthogonal and lie at right angles to each other. 

The unit vectors tangent to the loci are found using the following equation: 

Rt=-p==-X + -p=::ry 

JT+tt? Jl+ m^ ^ 

Where Rt = the unit vector tangent to the loci, and 

m = the slope of the loci, dyt/dxt at the emitter. 

The equations for the slope of the loci at an arbitrary point are developed in 


Appendix C and presented here: 


dyt _ [(xt - xa)^sin y - 2(xt - xa)(yt - ya)cos v|/ - (yt - ya)^sin y] 
dxt [(yt - ya)^cos vj/ - 2(xt - xa)(yt - ya)sin v}/ - (xt - xa)^cos v] 

Where xt, yt = the x and y coordinates of the emitter, 

xa, ya = the X and y coordinates of receiver A, 

xb, yb = the X and y coordinates of receiver B, and 
vf; = the bearing from receiver A to receiver B. 

The bearing from receiver A to receiver B is calculated from the following equation: 


V = arctan 


F yb-ya ' 

Lxb-xa. 




4. Ambiguity in TDOA Observations 


Ambiguity in the burst TDOA observations is present if another burst is received 
by one of the receivers before all of the receivers detect the first burst. In the burst 
TDOA problem each of the receivers is scaimed only once per rotation of the emitter. 
Ambiguity is possible if the distance separation between the receivers causes a burst fi-om 
the next scan of the emitter to be received by a closer receiver before a distant receiver 
detects the burst fi-om the first scan. This occurs if the time required for the signal to 
travel firom one receiver to the next is greater than the time required for the emitter to 
complete one scan. This is unlikely. In even the fastest scanning emitters the scan rate is 
unlikely to be greater than 10 revolutions per second. For an emitter with a scan rate of 
10 revolutions per second, the separation required to introduce ambiguity in the burst 
TDOA is 30,000 kilometers. 

B. SINGLE PULSE TIME DIFFERENCE OF ARRIVAL (TDOA) 

1. Problem Geometry 

If two sensors are illuminated at the same time by the beam of an emitter sending 
out pulsed signals, the pulse TDOA is the difference in the amount of time required for a 
single pulse to reach the two sensors. Since electromagnetic waves travel at 
approximately the speed of light, the TDOA is a function of the difference in the distance 
that the pulse must travel. The difference in distance the pulse must travel is calculated 
from the geometry of the problem. Figure 10 is a diagram of the geometry of the 
problem. 
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Figure 10. Single pulse TDOA geometry 


The time difference of arrival for a pulse at sensors A and B is calculated from 
equation (24). 


[(xt - xa)^ + (yt - ya)^ ] - [(xt - xb)^ + (yt - yb)^ ] 

TDOA =- *—^ - - — 

Where TDOA = the time difference of arrival, 

xt, yt = the X and y coordinates of the emitter, 

xa, ya = the x and y coordinates of sensor A, 

xb, yb = the X and y coordinates of sensor B, 
c = the speed of light. 


(24) 
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2. Loci of Constant TDOA 


Given a pulse TDOA observation and sensor geometry, the possible location of 
the emitter is not unique and a locus of possible emitter locations exists. Figure 11 is a 
plot of equation (24) for constant values of TDOA and the sensor configuration shown. 
The sensors are separated by 500 meters. For the TDOA observations shown, the 
location of the emitter lies anywhere along these hyperbolic loci. 



Figure 11. Loci of constant pulse TDOA, 
sensor separation 500 meters 


If more than two sensors are used it is possible to determine a unique estimate of 
the emitter's location. The scenario shown in Figure 12 uses two widely separated 
receiver platforms with two sensors attached to each of the receiver platforms. With this 
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configuration, the loci of constant TDOA intersect at a single point and a unique estimate 
of the emitter's location can be obtained from the pulse TDOA observations. 
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Figure 12. Loci of constant TDOA for two receiver platforms, 
sensor separation 500 meters 


3. Orthogonality of the Loci of Constant TDOA 

As previously developed in the burst TDOA problem the orthogonality of the 
TDOA observations can be related to the orthogonality of the loci of constant TDOA. 
The dot product of the unit vectors tangent to the loci give the cosine of the angle 
between the loci and a good indication of the orthogonality of the TDOA observations. 
The unit vectors tangent to the loci are calculated using equation (21). The equations for 
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the slope of the loci of constant TDOA are developed in Appendix C. The slope of the 
loci, dxt/dyt, is given by the following equation: 


dyt _ (yt - ya)^cos ly - (xt - xaXyt - ya)sin i}/ 
dxt (xt-xa)(yt-ya)cosv|/-(xt-xa)^sLn\|/ 

Where xt, yt = the x and y coordinates of the emitter, 

xa, ya = the X and y coordinates of sensor A, 

xb, yb = the X and y coordinates of sensor B, and 
v(/ = the bearing from sensor A to sensor B. 

The bearing frt>m sensor A to sensor B is calculated from the following: 




= arctan 


(yb-ya) 

(xb-xa) 


(25) 


(26) 


4. Ambiguity in Single Pulse TDOA Observations 

Ambiguity in pulse TDOA observations develops if a second pulse is detected by 
one of the sensors before both of the sensors detect the first pulse. This can occur if the 
amount of time required for a pulse to travel from one sensor to the next is longer than the 
PRI. Therefore, the separation of sensors determines the highest PRF that signals can 
have and still be detected without ambiguity. 

PRFm« = 

PRI = Rab (27) 

Where Rab = the distance between sensor A and sensor B, and 

c = the speed of light. 

For sensors separated by a distance of 500 meters the smallest PRI an emitter could have 
and be detectable without ambiguity is 1.66 microseconds. This corresponds to a PRF of 
approximately 600 kHz. 
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Estimates of the PRI of the emitter can be calculated from the time interval 


between pulses at each of the sensors. Therefore, the PRI can be estimated independently 
and used to test the TDOA observations for ambiguity. 
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rv. KALMAN FILTERING 


A. THE EXTENDED KALMAN FILTER 

There are several methods of linearizing a nonlinear problem to allow the Kalman 
filter to be used. The problem itself can be linearized by expanding the equations of state 
in a Taylor series around some norm. This approach is often used to derive equations of 
state in aircraft stability and control problems. The nonlinear equations of state are 
expanded and derived for pertubations around some equilibrium trim point. As long as 
the pertubations from the trim point are within acceptable limits this linear approximation 
will be acceptable and linear Kalman Filter theory can be q)plied. In the Extended 
Kalman Filter, the plant and full nonlinear equations of state are used, but the Kalman 
filter gain is calculated using a linear approximation of the model. The general equations 
of state are given by: 

x(k + 1) = f(x(k), w(k), k) (28) 

Where: x(k) = state of the system,. 

w(k)= the plant driving noise. 

The function f(«) can be a nonlinear function of the states, the noise, or k. 
The general equations for the observations are given by: 

z(k) = h(x(k),v(k),k) (29) 

Where: z(k) = the measurements of the system, and 

v(k) = the measurement noise of the system. 
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The function h(*) can be a nonlinear function of the states, the measurement noise, or k. 
The measurement noise and the plant driving noise are assumed to be zero mean white 

noise. 

The extended Kalman filter is given by equations (30) through (32). These equations 
utilize the nonlinear relationships to predict the states and observations of the the system. 
The smoothed estimates of the states are calculated using the Kalman gain calculated 
with a linearized version of the Kalman gain equation. 

x(k+llk) = f(x(klk),k) (30) 

z(k+llk) = h(x(k+llk),k) (31) 

x(k+llk+l) = x(k+llk)+G(k+l)(z(k+l)-z(k+llk)) (32) 

Where x(klk) = the current estimate of the states, 

x(k+llk) = the predicted estimate of the states, 
z(k + 1) = the observed measurements, 

G(k +1) = the Kalman gain, 

w(k), the plant driving noise, is assumed to be zero mean. 
v(k), the measurement noise, is assumed to be zero mean. 

Equation (30) is used to predict the next state of the system, given the current state of 
the system and the current input. Since the plant driving noise is assumed to have zero 
mean, it is assumed that the plant driving noise does not contribute the expected value of 
the next state. The next observation is predicted using equation (31). This equation uses 
the predicted state of the system. The measurement noise is assumed to have zero mean, 
and does not contribute to the the expected value of the observation. The smoothed 
estimates of the states of the system are found with equation (32). This equation uses the 
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Kalman gain ^^ch is calculated using equations (33) through (35). These equations use 
first order linear approximations of state prediction and observation equations. The 
covariance prediction equations and the Kalman gain equations are given as: 


P(k+llk) = 6P(k!k)6'f + Q 

(33) 

P(k +1 Ik +1) = [l - G(k + l)H]p(k +1 Ik) 

(34) 

G(k + 1) = P(k + 1 lk)HT[HP(k + 1 lk)HT + r]'' 

(35) 


Where P(klk) = the current estimate of the covariance of the estimation error, 

P(k+llk) = the predicted covariance, given the current value, 

P(k + 1 Ik + 1) = the predicted covariance given the current observation, 
Q = the covariance of the plant driving noise, 

R = the covariance of the measurement noise, 

4) = ^(x(k),u(k),w(k),k) 

x(k)=x(klk). u(k>=u(k), w(k)=0 

^ ^ gh(x(k),v(k),k) 

x(k>=x(klk). v(k)=0 


To initialize the Kalman filter, initial estimates for the states and the covariance 
matrix must be provided. The values chosen are generally: 

x(OlO) = E[x(0)], and (36) 

P(0l0) = E[{x(0)-x(0l0)}{x(0)-x(0l0)}T] = cov[x(0)] (37) 

The Kalman filter provides optimal performance, and guarantees convergence and 
stability, but the extended Kalman filter does not. In many applications the extended 
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Kalman filter provides accurate estimates of the states, but may not guarantee 
convergence and stability. Generally the convergence and stability are highly dependent 
on the values chosen for the initial estimates of the states, and covariance. In any 
application of the extended Kalman filter, the convergence and stability of the solutions 
must be thoroughly explored for as wide a range of initial conditions as is possible. 

B. ERROR ELLIPSOIDS 

The error ellipsoid is a means of visualizing the error in the estimate of the states of 
the Kalman filter. The uncertainty in the estimate is expressed in the error covarianc" P. 
The error covariance matrix P is the expected value of the covariance of the error in the 
states. 

P = E[(x-x)(x-x)'^] (38) 

Where x = the state estimate of the Kalman filter, 

X - the mean of the states of the Kalman filter 

If the noise entering the Kalman filter is Gaussian then the estimates of the states of 
the Kalman filter also have a Gaussian distribution. This is because a linear 
transformation of a set of Gaussian random variables will result in another set of 
Gaussian random variables. 
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The equation for the probability density function for a joint Gaussian random 


variable is; 


' (2n)’«id,. ; r ““ 


Where x = the vector of N random variables, 

X = the expected value of x, E[x], and 
P = the error covariance matrix of X. 

The error covariance matrix is given by the following: 


Pll 

P«2 

• • • Pin 

P21 

P22 

• . • P2N 

Pni 

PN2 

. . . Pnn 


(40) 


Where py = E[(xi - Xi)(Xj - xj)] 

A surface in N space can be found ^ere the probability of x will be a constant. This 
surface exists over the values of x for which the argument (x - x)^P“‘ (x - x) is a constant. 
The total probability that the values of x are lie on the surface defined is found by 


integrating the probability density function over the entire surface. These surfaces of 


constant probability are called error ellipsoids of constant probability. If the argument 
(x - x)^P~‘(x - x) is set equal to the constant c^ then the probabilities that the value of x 
will lie within the ellipsoid formed in N space is given in Table 1. 
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c 

1 

2 

3 

N 

1 

0.683 

0.955 

0.997 

2 

0.394 

0.865 

0.989 

3 

0.200 

0.739 

0.971 


Table 1. PROBABILITIES FOR ERROR ELLIPSOIDS 


In the Kalman filter application pursued here the estimates of the location of the 
emitter are joint Gaussian random variables. The mean of the steady state estimate of the 
location of the emitter is considered the true location of the emitter, and the difference 
between the estimate of the location of the emitter and the mean of the steady state 
estimates is the error. This error term is used to form the error covariance in 
equation (38), Redefining the error in the estimate of the location of the emitter, the 
argument of the joint probability density function is: 

c^=0I^P"'x) (41) 


, the error in the estimate of the emitter location, 

P = the error covariance matrix 
c = a constant. 

When the argument in equation (41) is expanded and the inverse of the P matrix is 
found, the following equation results: 


Pll Pl2 
Pl2 P 22 


, and 


Where 


X = 


xt 


?h-xt 

. 


^-yt 


c^ = 


1 


P 11 P 22 -P 12 


-[« yt ] 


' P 22 -P 12 1 

1- 

IX 

1 — 

1 

y 

— 1 

—j 


(42) 
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This further reduces to; 


P22^^ - 2 Puxtyt + Pnyt^ 

PnP22-PL 


= 


(43) 


To facilitate plotting this equation, the y estimate of the location of the emitter is solved 

for in terms of the x coordinates of the emitter. Equation (43) is rewritten in the form of a 
quadratic equation Ay^ + By + C = 0 wdiere 


A = 


II 


PiiP22-PT2 


B = 


-2P,2Xt 


C- 


PnP22-P|2 P11P22-P12 


The quadratic formula is used to calculate the value of ytin terms of xt and the elements 
of the error covariance matrix. The following equation results: 


yt = 


Pl 2 Xt^ |(P|lP 22 -Pj 2 )^ _2 


11 


-(Piic^-jh^) 


(44) 


11 


The maximum and minimum values of the error xtare found from the locations where 
xt^ = Piic^. If xt^ > Pnc^ a real solution to equation (44) will not exist. 
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V. SIMULATION RESULTS 
A. THE BURST TDOA KALMAN FILTER 


1. Extended Kalman Filter Equations 

An extended Kalman filter algorithm is developed to estimate of the emitter 
location from burst TDOA observations. Simulated data is used to test the accuracy and 
stability of the algorithm. The TDOA measurement noise is assumed to be white 

gaussian. The following assumptions are made to simplify the problem: 

1. The emitter scans in azimuth at a constant rate. 

2. An estimate of the scan rate of the emitter is known a priori. 

3. Some method is used to calculate the TDOA of the radar signal between the two 
receivers, and the error on this TDOA observation is modeled as additive white 
Gaussian noise with variance R. 

4. The emitter is stationary. 

5. The X and y coordinates of the emitter and receivers are measured with respect to 
a fixed coordinate reference. 

6. The exact locations of the receivers are known for all points in time. 

These assumptions simplify the problem, but do not reduce its usefulness, because the 

algorithm can easily be modified to account for these assumptions. 

The estimated state and observations are calculated from the following: 

x(k+l) = x(k)= (45) 


z.b(k) = 


J_ (xt - xa)(yt - yb) - (yt - ya)(xt - xb) 

[[(xt - xa)^ + (yt - ya)^][(xt - xb)^ + (yt - yb)^]] 


+ V|*(k) 


(46) 
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+ VkOc) (47) 


^ J_ (xt -xaXyt-yc)-(yt- ya)(xt - xc) 

[[(xt - xa)^ + (yt - ya)^][(xt - xc)^ + (yt - yc)^]] 

i 1 (xt-xb)(yt-yc)-(yt-ybXxt-xc) , 

Zbc(K) = — -;-;;-;-—— + Vbc(k) (48) 

[[(xt - xhy + (yt - yb)^][(xt - xc)^ + (yt - yc)^]] J 

Where xt,yt = the estimate of the emitter's x and y coordinates, 

SR = the scan rate of the emitter, 

v„ = the measurement noise present in the TDOA observation 
between two receivers, 

xa, ya = the X and y coordinates of receiver A, 

xb, yb = the X and y coordinates of receiver B, 

xc, yc = the X and y coordinates of receiver C, 

xt, yt, xa, ya, xb, yb, xc, and yc are all functions of k. 

Since the emitter is stationary, the state transition matrix in equation (45) is the 
identity matrix. Therefore, the predicted estimate of the emitter's location is the same as 
the current estimate of the location. Since the coordinates of receivers A, B, and C are 
known at all points in time and the scan rate is known a priori, the predicted TDOA 
observations are only functions of the predicted states. 

The Kalman filter is implemented as three separate Kalman filters linked 
together. Each receiver processes the received observations separately. They calculate 
Kalman gains and an error covariance matrix from its observations, and then use these to 
calculate an estimate of the location of the emitter. This estimate is then passed to the 
next receiver. This algorithm is used because it facilitates implementation in hardware. 
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The equations necessary for the Extended Kalman Filter implementation are: 


P(k+llk) = P(klk) + Q (49) 

P(k +1 Ik +1) = [l - G«(k + l)H«]p(k +1 Ik) (50) 

Gxx(k +1) = P(k +1 lk)HT [HxxP(k +1 \k)Hl, + V** ]"' (51) 

Where P(k+l|k) = the predicted error covariance, 

P(k|k) = the current error covariance, 

Q = the covariance of the state disturbances, 

G„(k+1) = the Kalman Gain calculated for the TDOA observation 
between tv:o of the receivers, 

V„ = the covariance of the measurement noise for the TDOA 
observation between two of the receivers, 

„ ^ r ai„(k) 3ho(k) I 
"w 1_ 3yt J 

The algebraic equations that make up the partial derivatives in H„ are extensive and are 
omitted here fo»- cb The derivation of is included in Appendix D. 

The flow chart for the Kalman filtering algorithm is shown in Figure 13. This 
flowchart provides the basic logic flow of the Kalman filtering algorithm. In an actual 
implementation, TDOA observations are available, and only the portion contained within 
the dotted outline is required. The additional portion of the algorithm generates the 
observations used to test the algorithm. 
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2. Burst TDOA Simulations 


The following scenario is used to test the burst TDOA Kalman filter algorithm. 
Three receivers are used and all three of the available TDOA observations are used to 
estimate the location of the emitter. The initial locations of the receivers and their 
velocities and bearings are shown in Figure 14. 


> 

< 

\ 

Receiver A 
(0.0) km 

^ V = 30 m/sec 

Bearing 90 deg. 

w V 

Receiver C 
(-20.-20) km 

V = 42 m/sec 

Bearing -135 deg. 

/ 

\ 

Receiver B 
(20.-2) Km 

V * 42 m/sec 

Bearing -45 deg 


Figure 14. Initial location and velocity of receivers 


The emitter locations are shown in Figure 15. 
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Although the aigorithm has the capability to account for different signal strengths 
and resulting error variances, all of the receivers are considered identical with identical 
signal to noise ratios. An 8 MHz sampling rate is assumed. The algorithms developed in 
Appendix A are used to calculate the variance of the error present in the TDOA 
observations for the signal to noise ratio and the 'ompling rates chosen. In the tests of 
Kalman filtering algorithm, this error is modeled as zero mean white Gaussian noise, that 


43 





is added to the true TDOA observations. The variance of the TDOA observation noise is 


5.660 X 10 ’ for all of the scenarios. 

For all of the scenarios the a priori estimate of the emitter location is chosen as 
(0,10) kilometers. The filter parameters Q, the covariance of the state excitation, and Po, 
the initial error covariance, are varied until the filter converges well. All of the scenarios 
are run for twenty five TDOA observations. Given the scan rate of the emitter, the 
real-time length of the simulation is twenty five seconds. 

Three plots are presented for each scenario. The first plot is an X-Y plot of the 
estimates of the coordinates of the emitter. This plot demonstrates the track of the final 
estimate of the emitter location after all three TDOA observations are processed for each 
time step. In this plot the loci of constant TDOA that correspond to the final steady state 
TDOA observations are also plotted. Some of the loci do not intersect exactly at the 
emitter. The equations used to plot the loci assume that the distance to the emitter is 
much greater than the separation of the receivers, and in some of the scenarios this 
assumption is not valid. The loci visually illustrate the relationship between the TDOA 
observations tmd give a visual indication of the orthogonality of the observations. The 
3o error ellipsoids are also plotted for the 1st, 8th, 15th, and 22nd estimates of the emitter 
location. These ellipsoids provide a visual indication of the accuracy of the location 
estimate. 


44 



The second plot is the estimate of the x and y coordinates of the emitter vs. the 
number of TDOA observations. This plot gives a visual indication of how rapidly the 
estimates of the emitter location converge on the true coordinates of the emitter. 

The third plot is an X-Y plot that shows the true emitter coordinates, the loci of 
constant TDOA, and the estimates of the emitter coordinates in the vicinity of the mean 
steady state estimate of the emitter coordinates. The 3o error ellipsoid is plotted for the 
steady state estimate, and the length of the major and minor axis of the error ellipsoid 
give an indication of the maximum error in the estimate. 
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a. Scenario U1 


In the first scenario the performance of the burst TDOA filter is tested for an 
emitter at a range of SOO kilometers from the origin and a bearing of 90 degrees from the 
X coordinate axis. The initial error covariance matrix, Po, and the variance of the state 
excitation, Q, are varied until the filter converges to a steady state value within ten 
observations. The plots presented in figures 16 through 18 are representative of the 
results obtained. The filter converges very well, and as shown in figures 16 and 18 the 
estimate of the location of the emitter converges to the true location of the emitter. The 


values that gave the best convergence were: 


Po = 


1.0 0 
0 1.0 


4.0 0 
0 4.0 


X 10"m^ 


xl0“m^and Q = 

As shown in Figure 16, the 3a error ellipsoids grow larger as the estimate of 
the location of the emitter grows closer to the actual location. Because the distance to the 
emitter is large in relation to the separation of the receivers, the TDOA observations are 
very close to each other. The estimated location of the emitter does not approach the true 
location until the error covariance matrix and subsequently the Kalman gains grow large 
enough to amplify the very small differences in the TDOA and drive the states of the 
filter. A high value for Q is required to drive the error covariance matrix high enough. If 
a smaller value of Q is used, the filter reached a steady state value far from the true 
location of the emitter. 
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Figure 16. Scenario #1 Estimates of the location of the emitter 
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Figure 18. Scenario #1 Close-up of steady state estimate of emitter location 
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b. Scenario if2 


In the second scenario, the performance of the burst TDOA filter is tested at 
an intermediate range. The emitter location is chosen at a range of 150 kilometers from 
the origin and at a bearing of 30 degrees from the x coordinate axis. The coordinates of 
the emitter are (130,75) kilometers. The initial error covariance matrix, Po, and the 
variance of the state excitation term, Q, are varied until the filter converges to a the 
correct location within approximately ten observations. The plots presented in figures 19 
through 21 are representative of the results obtained. The values that give the best results 
are: 


Po 


1.0 0 
0 1.0 


5.0 0 
0 5.0 


X lO^m^ 


xl0^°m^and Q = 

The filter performed very well. The orthogonality of the TDOA observations 
in the vicinity of the emitter is high, and the estimate of the emitter location converged to 
the proper location quickly and smoothly. The final steady state estimate is stable and no 
numerical problems are present. The 3a error ellipsoids are not visible in figure 19 
because of the scale of the figure. 
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Figure 19. Scenario #2 Estimates of emitter location 
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Figure 20. Scenario #2 Estimates of emitter coordinates vs observations 



lop -------p---^ 

i ; y 

1 Major Anas 

3.1 km . .y 

77| Minor Axis* 

1.2 km.' // 

, 




1 

1 

/ " ^ j ''' ' 

7»; Enww _ 


i ^ 


74,,--' /V 


1 

, >9 error ciiipaoid 

73 

\ 


loci of mnini TBOA j 


X coordinate (kilometers) 



Figure 21. Close-up of the steady state estimate of emitter location 











c. Scenario #5 


In the third scenario the performance of the burst TDOA filter is tested for an 
emitter at close range. The emitter location is chosen at a range of 30 kilom^ers fi-om the 
origin and at a bearing of 60 degrees fi-om the x coordinate axis. The initial error 
covariance matrix, Po, and the variance of the state excitation term, Q, are varied imtil the 
filter converged to the correct location of the emitter within approximately ten 
observations. The plots presented in figures 22 through 24 are representative of the 


results obtained. The values that give the best results are; 


1.0 0 
0 1.0 


X 10‘* m^ and Q = 


40 0 
0 40 


Although the loci of constant TDOA plotted in figure 22 do not intersect at 


the emitter, they indicate that the orthogonality of the TDOA observations is high and 


accurate results are expected. Figure 23 shows that the estimates of the emitter location 


converged to the correct location of the emitter rapidly and smoothly. Because of the 


large scale of the plot, the 3a error ellipsoids are not visible in figure 22. 
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Figure 22. Scenario #3 Estimates of emitter location 







Figure 23. Scenario #3 Estimates of emitter coordinates vs observations 



Figure 24. Close-up of the steady state estimate of emitter location 
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3. The Burst TDOA Filter Results 


The simulations presented here in no way test all the possible implementations of 
this filtering algorithm. Further evaluation of the algorithm is necessary to fully evaluate 
its performance and capabilities given the wide variety of possible emitter locations, 
receiver configurations, and emitter characteristics. From the an<ilysis conducted some 
general observations can be drawn. 

For the burst TDOA Kalman filter the orthogonality among the TDOA 
observations is generally very good. For the three randomly chosen receiver locations, 
the orthogonality of the TDOA observations are not heavily dependent on the location of 
the emitter. As the locations of the receivers remained constant and the bearing to the 
emitter is varied throughout a 90 degree arc, and there does not appear to be any 
particular bearings where the orthogonality of all three observations would decrease to 
the point where the performance of the filter degrades. For all bearings and ranges 
examined at least two of the TDOA observations had good orthogonality. 

The accuracy of the filter is heavily dependent on the values chosen for Q, the 
covariance of the state excitation. If this value is too low, the filter reaches a steady state 
value, but the steady state value is not the correct location of the emitter. As Q is 
increased, the steady state estimate of the location of the emitter reaches a stable point 
around the true emitter location. The magnitude of Q directly affects the magnitude of 
the error covariance P, and as Q is increased the error covariance and the 3a error 
ellipsoids for the steady state estimate increase in size. 


55 





The filter performs well for the receiver configurations and emitter locations 
tested. Even at the maximum range, where the receivers are spaced closely in comparison 
to the distance from the emitter, the error present in the final estimate of the location of 
the emitter is reasonable. For scenario #1 the length of the 3a error ellipsoid, which 
represents effectively the maximum error in the estimate given the noise present in the 
observations, is roughly 30 kilometers in bearing direction and 10 kilometers in range. 
These represent errors of about 3 degrees in bearing and 2 percent in range. In scenarios 
#2 and #3, the distance to the emitter decreases relative to the separation of the receivers 
and the maximum error decreases to about 0.5 and 0.20 degrees in bearing, and 2 percent, 
and 0.3 percent in range respectively. 

These results are accomplished filtering the TDOA observations one at a time. 
This filtering technique tends to skew the error ellipsoids to align them vrith the loci of 
constant TDOA that corresponded to the last TDOA observation. This places the major 
axis of the error ellipsoid along this loci and tends to increase the error along that loci. If 
all three of the TDOA observations are processed simultaneously as the filter reaches a 
steady state solution, the error covariance and thus the error ellipsoids and maximum 
error present in the estimate can be decreased further. 
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B. SINGLE PULSE TDOA KALMAN FILTER 


1. Extended Kalman Filter Equations 

An extended Kalman filter algorithm is developed to estimate of the emitter 
location fix>m pulse TDOA observations. Simulated data is used to test the accuracy and 

stability of the algorithm. The following assumptions are made to simplify the problem: 

1. Two receiver platforms are used, and each receiver platform is equipped with two 
sensors to detect the pulse TOA. 

2. The error present in the pulse TDOA observation is modeled as zero mean 
additive white Gaussian noise with variance R. 

3. The emitter is stationary. 

4. The X and y coordinates of the emitter and receivers are measured with respect to 
a fixed coordinate system, and the coordinates of the sensors are known exactly. 

5. The PRI is long enough that ambiguity is not present in the pulse TDOA 
observations. 

These assumptions simplify the problem, but do not reduce its usefulness, because the 
algorithm developed can easily be modified to account for these assumptions. 

The estimated state and predicted observation are calculated fi'om the following: 

x(k+l) = x(k)= (52) 


z(k+llk) = 

The Kalman filter equations are: 

x(k+llk+l) = x(k+llk) + G(k+l)[z(k+l)-z(k+llk)] (54) 

P(k+llk) = P(klk) + Q (55) 


(xt(k) - xa(k))^ + (yt(k) - ya(k))^ 1 - [ (xt(k) - xb(k))^ + (yt(k) - yb(k))^ 
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(56) 


P(k + 1 Ik 1) = [l - G(k + l)H(k)]p(k +1 Ik) 

G(k + 1) = P(k + 1 lk)H(k)'r[H(k)P(k + 1 lk)H(k)T + r]'‘ (57) 

Where xt(k) ^(k = the estimate of the x and y coordinates of the emitter, 

z(k+l) *= the TDOA observations for receiver platforms 1 and 2, 
z(k +1 Ik) = the predicted TDOA based upon the estimates of emitter 
location. 

xa, ya = the X and y coordinates of the .sensor A of the receiver platform, 

xb, yb = the X and y coordinates of the sensor B of the receiver platform, 
G(k+1) = the Kalman gains, 

P(k|k) = the error covariance, 

H(k) = the gradient of the observation equation, 

Q = the covariance of the plant noise, 

R = the covariance of the measurement noise, 

c = the speed of light, 

xa, xb, ya, yb are all functions of k. 

The gradient of the observation equation H(k) is calculated from the following equations. 



TT/n.x r lf(xt(khxa) (xt(k>-xb)'l lf(yt(k)-ya) (yt(k>-yb)'\ 

cV Ra(k) Rb(k) J Raflt) Rbflc) J 

(58) 

Where, 

Ra(k) = [(5rt(k) - xa)^ + (^(k) -- ya)^ J ^ , and 

(59) 


Rb(k) = [(xt(k) - xb)^ + (>h(k) - yb)^ ] 

(60) 


The flowchart for the pulse TDOA Kalman filter algorithm is identical to the 
flowchart used for the burst TDOA Kalman filter shown in Figure 13. The simulation 
calculates the pulse TDOA observation using the known locations of the emitter and 
sensors. Zero mean white Gaussian noise is ar* Jed to the calculated observations and the 
filter estimates the location of the emitter from these noisy observations. The variance of 
the measurement noise is calculated from the emitter parameters and an assumed signal to 
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noise ratio. The filtering algorithm is tested for a number of emitter locations. In each 
scenario the initial error covariance matrix and the variance of the plant excitation noise 
are varied until the filter converges to the emitter's true coordinates 

2. Pulse TDOA Kalman Filter Simulations 

Three scenarios are used to test the performance of the pulse TDOA Kalman 
filter. The receiver locations remain the same for all of the scenarios and are shown in 
Figure 25. The sensors are stationary for all of the scenarios. The sensors are separated 
by 500 meters, and the receiver platform separation is 20 kilometers. 
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Figure 25. Receiver platform and sensor configuration 
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The emitter locations chosen to test the pulse TDOA algorithm are shown in Figure 26. 



Figure 26. Emitter locations for the single pulse TDOA Kalman filter 


The emitter has the following characteristics. These characteristics are typical for long 


range search radar. 

Beam width 
Scan Rate 
PRF 

Pulsewidth 

Peak Signal to Noise Ratio 


2.0 deg. 

360 deg/sec. 
2000 Hz 
1.0 microsecond 
15dB 


All of the sensors are considered identical with identical peak signal to noise 
ratios. An 8 MHz sampling rate is assumed. The algorithms developed in Appendix A 
are used to calculate the variance of the error present in the TDOA observations for the 
sampling rate and signal to noise ratio chosen. In the tests of the Kalman filtering 
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algorithm, this error is modeled as zero m jan white Gaussian noise, and is added to the 
true TDOA observations. The variance of the TDOA observation noise is 1.694 x 10 '* 
for all of the scenarios. 

For all of the scenarios the a priori estimate of the emitter's location is 
(10,5) kilometers. This a priori estimate provided good convergence for a wide range of 
emitter locations. The variance of the plant excitation noise Q, and the initial error 
covariance are varied until the filter converges to the actual coordinates of the emitter. 
All three of the scenarios are run for sixty TDOA observations. For the emitter 
characteristics chosen, this represents a real run time of 30 milliseconds. 

Three plots are presented for each of the scenarios examined. The first plot is an 
X-Y plot of the estimates of the coordinates of the emitter. This plot demonstrates the 
track of the final estimate of the emitter location after both TDOA observatioas are 
processed at each time step. The loci of constant TDOA that correspond to the final 
steady state TDOA observations are also plotted. These loci illustrate the relationship 
between the TDOA observauons and give a visual indication of the orthogonality of the 
observations. Also, the 3a error ellipsoids are plotted for the 1st, 20th, 39th, and 58th 
estimates of the emitter location. These ellipsoids provide a visual indication of the 
accuracy of the location estimate. 

The second plot is the estimate of the x and y coordinates of the emitter vs. the 
number of TDOA observations. This plot gives a visual indication of how rapidly the 
states of the filter converge to the true coordinates of the emitter. 
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The third plot is an X-Y plot of true emitter coordinates, the loci of constant 
TDOA, and the estimates of the emitter coordinates in the vicinity of the mean steady 
state estimate of the emitter coordinates. The 3a error ellipsoid is plotted for the steady 
state estimate. The length of the major and minor axis of the error ellipsoid give an 
indication of the maximum error in the estimate. 
a. Scenario U1 

In the first scenario the performance of the pulse TDOA filter is tested for an 
emitter located at a range of 500 kilometers fiom the origin at a bearing of 90 oegrees 
fi-om the X coordinate axis. The initial error covariance matrix, Po, and the variance of 
the state excitation, Q, are varied until the filter converges to a steady state value. The 
plots presented in Figures 27 through 29 are representative of the results. In each 
simulation, a total of sixty observations are filtered. The values that give the best results 
are: 


Po = 


5.0 0 
0 5.0 


xlO^^m^and Q = 


1.0 

0 


0 

1.0 


X 10* m^ 
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Figure 27. Scenario#! Estimates of emitter location 





















The filter performs well, but the error in the steady state estimate is high, and 
the estimate of the y coordinate of the emitter, shown in Figure 28, converges very 
slowly. Because the distance to the emitter is large in relation to the separation of the 
receivers, the TDOA observations are very close to each other. The estimated location of 
the emitter does not approach the true location of the emitter until the error covariance 
matrix and subsequently the Kalman gains grow large enough to amplify the very small 
differences in the TDOA to drive the states of the filter near the true location. A large 
value for Q is required to drive the error covariance matrix high enough. If a smaller 
value of Q is used, the filter reaches a steady state far fi-om the true location of the 
emitter. Because of the magnitude of the error covariance matrix, the 3o error ellipsoids 
grow large. The major axis of the steady state error ellipsoid is nearly 400 kilometers. 
Considering that the distance to the emitter is only 500 kilometers, a location estimate 
with that much possible error has limited usefulness. 

In a second attempt, the separation of the receiver platforms is increased to 
300 kilometers to determine if greater separation of the receiver platforms can increase 
the orthogonality of the TDOA observations and improve the performance of the filter. 
The separation of the sensors is maintained at 500 meters. The results of the simulation 
are shown in Figures 30 through 32. The value chosen for Q in this simulation is 
x 10*’. The value chosen for Po is the same as the previous simulation. 


' 2.0 0 
0 2.0 
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Figure 31. Scenario #1 Estimates of emitter coordinates vs. observations for 

300 km receiver separation 



Figure 32. Scenario #1 Close-up of steady state estimate of emitter location for 300 km 

receiver separation 
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As shown in Figure 31 the estimates of the emitter coordinates approached the 
true coordinates of the emitter more quickly and directly than the previous simulation. 

This is largely due to the greater oithogonahty of the TDOA observations. The steady 
state 3a error eUipsoid plotted in Figure 32 is considerably smaller than in the previous 
simulation and the major axis of the error eUipsoid is only S1 kilometers. 
b. Scenario U2 

In the second scenario the performance of the burst TDOA filter is examined 
with the emitter located at an intermediate range. An emitter is located at a range of 1 SO 
kilometers from the origin and a bearing of 30 degrees fi^om the x coordinate axis. The x 
and y coordinates of the emitter are (130, 75) kilometers. As in the first simulation, the 
initial error covariance matrix, Po, and the variance of the state excitation term, Q, are 
varied until the filter converges to the true coordinates of the emitter. The values that give 
the best convergence are: 


Po = 


1.0 0 
0 1.0 


xlO m , and Q = 


1.0 0 
0 1.0 


xl0’m^ 


The plots presented in Figures 33 through 35 are representative of the results 


obtained. 
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Figure 33. Scenario #2 Estimates of emitter location 
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Figure 34. Scenario #2 Estimates of emitter coordinates vs. observations 



Figure 35. Close-up of the steady state estimate of emitter location. 
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The filter performs well, but the error in the steady state estimate is high, and 
the estimates of the x and y coordinates of the emitter, shown in Figure 34, converge very 
slowly. The distance fi’om the receiver platforms to the emitter is much less than in 
scenario #1, but the orthogonality of the TDOA observations is low because of the 
location of the emitter in relation to the receivers. As shown in figures 33 and 35, the loci 
of constant TDOA are nearly parallel in the vicinity of the emitter. Although the filter 
converges to a steady state value close to the actual emitter location, the error covariance 
is very high. The major axis of the 3o error ellipsoid for the steady state estimate is 
approximately 140 kilometers. Considering that the distance to the emitter is only 150 
kilometers, a location estimate with an error ellipsoid this large has limited usefulness. 

As shown in the first scenario, greater separation of the receiver platforms can 
increase the orthogonality of the TDOA observations and improve the performance of the 
filter. When the separation of the receiver platforms is increased to 150 kilometers and 
the separation of the sensors maintained at 500 meters, the performance improves 
dramatically. The results of the simulation are shown in figures 36 through 38. The final 

r 5.0 0 1 , 

value chosen for Q is q 5 q ^ value chosen for Po is the same as the 

previous simulation. 
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Figure 36. Scenario #2 Estimates of the location of the emitter with 
150 km receiver separation 
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Figure 38. Close-up of the steady state estimate of emitter location with 150 km receiver 

separation 
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As shown in Figure 37, the estimates of the emitter coordinates approach the 
true values more quickly and directly than the previous simulation. This is largely due to 


the greater orthogonality of the TDOA observations. The 3o error ellipsoid plotted in 
Figure 38 is considerably smaller than in the previous simulation The major axis of the 
error ellipsoid is only 16.7 kilometers, 
c Scenario U3 

In the third scenario the performance of the pulse TDOA filter is tested for an 
emitter located at close range in comparison to the separation of the receivers. The 
emitter location is located at a range of 30 kilometers fi'om the origin and at a bearing of 
60 degi ces from the x coordinate axis. The x and y coordinates of the emitter are (15, 26) 
kilometers As in the previous simulations, the initial error covariance matrix and the 
vanance of the state excitation term, are varied until the filter converges to the true 
coordinates of the emitter. The values that give the best convergence are: 


Po = 

1 

o 

o 

1_ 

xl0'*m^and Q = 

1- 

o 

o 


} 

o 

o 

_1 


L 0 40 J 


The plots presented in Figures 39 through 41 are representative of the results 

obtained. 

As shown in Figure 40 the estimate of the emitter location converges to the 
true coordinates of the emitter quickly. Increasing the value of Q improves the 
convergence of the filter, but also increases the error present in the final estimate and 
increases the size of the error ellipsoid For the values of Q chosen the major axis of the 
3o error ellipsoid is only 1.4 kilometers 
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Figure 39. Scenario #3 Estimates of the emitter location. 
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Figure 40. Scenario #3 Estimates of emitter coordinates vs. observations 



Figure 41. Close-up of the steady state estimate of the emitter location 
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3. Pulse TDOA Filter Results 


The simulations presented here do not test all the possible implementations of 
this filtering algorithm. Further evaluation of the algorithm is necessary to fully evaluate 
its performance and capabilities given the infinite variety of possible emitter locations, 
receiver configurations, and emitter characteristics. From the analysis conducted some 
general observations can be drawn. 

The pulse TDOA filter accuracy and effectiveness is heavily dependent on the 
orthogonality of the TDOA observations. As the orthogonality of the TDOA 
observations decreased, the error present in the estimate increased. When the receiver 
platforms are closely spaced relative to the distance to the emitter, the loci are nearly 
parallel to each other and the TDOA observations have low orthogonality. To obtain 
good filtering results the spacing between the receiver platforms has to be increased. 

The receiver locations chosen for the single pulse TDOA filtering problem 
yielded areas where the orthogonality of the TDOA observations is low. In these areas, 
the ability of the algorithm to estimate accurately the location of the emitter is limited. 
These blind areas can be avoided with different sensor configurations or more receivers. 

The accuracy of the algorithm filter is also heavily dependent on choosing the 
correct values for the covariance of the state excitation, Q. If this value is too low, the 
filter reaches steady state, but the steady state estimate of the location is far from the true 
location of the emitter. As Q is increased, the algorithm reaches a steady state around the 
emitter's true location, but the error covariance and the error present in the final estimate 


77 



rises much higher. Consequently the size of the 3a error ellipsoids and the maximum 
error present in the final estimates are much larger. 

The final approaches used in each of the scenarios yield good results. In the 
second approach used in scenario #1, the accuracy obtained from the filtering is 
approximately 3 degrees in bearing and about 10 percent range. For scenario #2 the error 
is approximately 5 degrees in bearing and 5 percent in range, and in the third scenario, the 
error is approximately 2 degrees in bearing and 3 percent in range. 

The pulse TDOA algorithm filters the TDOA observations one at a time. This 
filtering technique tends to skew the error ellipsoids and align them with the loci of 
constant TDOA that corresponded to the last TDOA observation. This places the major 
axis of the error ellipsoid along this loci and tends to increase the error in that direction. 

If both of the TDOA observations are processed simultaneously as the filter reaches 
steady state, the error covariance and the error ellipsoids can be decreased. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 

Both of the Kalman filtering algorithms considered here perform well. Given the 
proper parameters, both algorithms accurately estimate the location of the emitter to a 
reasonable degree of accuracy. The orthogonality of the TDOA observations is the 
largest factor that affects the accuracy of the final estimates. 

Overall, the geometry of the burst TDOA filtering problem yields TDOA 
observations with better orthogonality. Even when the distance to the emitter is much 
larger than the distance between the receivers, the orthogonality of the TDOA 
observations are good. In the burst TDOA filtering problem, the orthogonality of the 
observations is not as sensitive to the location of the emitter in relation to the receivers. 
The receiver configuration provides good orthogonality for a vride range of emitter 
locations. To obtain good orthogonality in the pulse TDOA filtering problem, the 
receivers require wide separation. 

The major disadvantage of the burst TDOA filter is the slow rate at which 
observ'ations are obtained. Observations are obtained at the scan rate of the emitter, and 
for the scenarios examined, this means that only one observation is made per second. For 
the simulations pursued, the total run time is 25 seconds, but reasonably accurate 
estimates of the location are available within about 10 observations. The number of 
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observations required to reach steady state could be reduced further with better a priori 
estimates. 

The most attractive feature of the pulse TDOA filter is the rate at which observations 
are obtained. Observations are available at the rate of the PRF of the emitter, and for the 
scenarios pursued here, the pulse TDOA filter could obtain 2000 observations for every 
observation obtained by the burst filter. 

It may be possible to combine the two filtering algorithms and utilize the best 
features of both. As stated previously, with only two receivers the location of the emitter 
can not be determined uniquely from the burst TDOA observations. However, if the 
receivers are equipped with multiple sensors, as in the pulse TDOA filtering problem, 
bearing estimates could be obtained and the location of the emitter along the burst TDOA 
locus could be uniquely determined. 

B. RECOMMENDATIONS 

The analysis and development conducted here is limited in its scope. Further testing 
and evaluation are required to more accurately evaluate the algorithms and models 

presented. Specifically, the following are recommended 

1. Further evaluation and testing of the burst and pulse TDOA algorithms be done to 
assess their performance for a wider variety of receiver and emitter configurations 
and scenarios. 

2. Perform a detailed analysis of the effect of a priori estimates on the convergence 
and accuracy of both filters. 

3. Explore the possibility of using the orthogonality of the TDOA observations to 
pick the best observations to filter and its effect on the performance of the filter. 

4. Obtain more detailed models of the error present in the burst and pulse TDOA 
observations, and determine their effect on the performance of the filter. 

5. Explore the possibility of combining the two filtering approaches to take 
advantage of the best features of both. 
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APPENDIX A. TOA PROBABILITY DENSITY 
A. PULSE TOA PROBABILITY DENSITY 

1. Derivation of Probability Density 

The pulses received and detected by the sensors are assumed to have the 
envelope shown in Figure A-1. 



Figure A-1. Pulse envelope and sampled pulse 


The sampling interval is assumed much smaller than the pulse width so the sampled 
pulse is a reasonable representation of the original pulse. The higher the sampling 
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probability density function of the pulse TO A, each sample of the pulse is considered a 
random variables as shown in equation (A-1). 

2(k) = Zo(k) + v(k) (A-1) 

Where z(k) = the random variable representing the amplitude of the kth sample, 

Zg(k) = the deterministic variable that is the true amplitude of the kth 
sample, and 

v(k) = the white Gaussian random variable that represents the noise 
present in the amplitude of the kth sample. The noise has zero 
mean and a variance of V, i.e.. E[v(k)^] = V. 

The mean and variance of z(k) are; 

Pz(k) = E[zo(k) + v(k)] = E[Zo(k)] + E[v(k)] = Zo(k) (A-2) 

= El (z(k) - P 2 (k))(z(k) - pz(k)) ] (A-3) 

Ojoc) = E[ (zo(k) + v(k) - Zo(k))(z<,(k) + v(k) - Zo(k)) ] = E[v(k)2] = V (A-4) 
Each sample is uncorrelated. 

E[z(k)z(k + t)] = E[ZoOt)Zo(k + x)] = 0 (A-5) 

Since they are Gaussian random variables, they are independent. 
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Equation (A-6) is used to calculate the time of arrival for the centroid of the pulse. 


Y I t(k)2(k) 

TOA = § = !^=i5- (A-6) 

Where TOA = the time of arrival for the pulse centroid, 

X = numerator of the pulse centroid TOA equation, 

Y = denominator of the pulse centroid TOA equation, 
t(k) = the time at which the kth sample of the pulse is taken, 
z(k) = the amplitude of the kth sample of the pulse, and 
N = the number of samples in the pulse. 

The time base from which the time of the samples are measured is arbitrary, but is 
common for all TOA measurements. This value is assumed deterministic and is therefore 
known exactly for each sample. 

The numerator and denominator of equation (A>6) are themselves Gaussian random 
variables. As shown in reference [1], if W is a linear combination of independent 
Gaussian distributed random variables. 


W = ai X(l) + a2X(2) + a3X(3) +.a„X(n) (A-7) 

then W i.s a Gaussian distributed random variable with a mean and variance given by; 

HW = ai|iX(l) + a2tlX(2) + a3tlx(3) +.anttxcn) (A-8) 

I 

I tJw = ai^^l)+ ^20^2) ®3tT^3) + • •antJ^n) (A-9) 


Based upon these relations the mean and variance of the numerator and denominator of 
equation (A-6) are given by; 
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(A-10) 


\ix = I t(k)z„(k) 

k»l 


oJ=It(k)'V (A-ll) 

k«i 

^Y=2^Zo(k) (A-12) 


o^= I V = NV 

k=l 


(A-13) 


Equation (A-6), the estimate of the TOA of the pulse centroid, is the ratio of two 
Gaussian distributed random variables. The joint probability density function for two 
Gaussian distributed random variables x and y is; 


exp 


K\y) = 


2(l-pi) 




2jto*ay(l -p^) 


2\l/2 


Where p, = the mean of x and y, 

o, Oy = the standard deviation of x and y, 
p = the correlation coefficient which is given by; 


(A-14) 


P = 


cov[xy] E[(x- p,)(y- py) 

OxOy OxOy 


(A-15) 


For equation (A-6), the covariance and correlation coefficient are calculated from the 
following: 


cov[X, Y] = E[(S t(k)z(k) - i t(k)zo(k))(i z(k) - S Zo(k))l (A-16) 

k«l k=l k=l k=l 


cov[X, Y] = E[(|^ t(k)v(k))(^I v(k))] 


(A-17) 
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Since the noise term v(k) is considered white noise with zero mean and variance V, the 
cross terms in the multiphcation in equation (A-17) will be zero and the resulting 
covariance is; 

cov[X, Y] = t(k)v(k)2 = t(k)V (A-18) 

Thus the correlation coefficient is: 

VEt(k) Et(k) 

n =- — -= —- rA-lQi 

The correlation coefBcient does not depend upon the variance of the noise V only on the 
number of samples and the sampling interval. 

From reference [2], the probability density function, f^(z), for the raho of two random 

y 

variables, Z = ^, is found from the following equation: 

f^(z)= I |y|f(zy,y)dy (A-20) 

-«0 

Where f(zy,y) is the joint probability density function f(x,y) with the variables zy 
substituted for x. 
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The probability density flmction for the pulse centroid TOA is: 

i - 2.a.a,(.-pv^ -^ 


(A-21) 


Where z = the pulse centroid TOA, 

y = the denominator of the pulse centroid TOA equation, 

= the mean of the numerator and denominator of the centroid TOA 
equation, calculated from equations (A-10) and (A-12), 
a/ Oy = the variance of the numerator and denominator of the centroid 

TOA equation, calculated from equations (A-11) and (A-13), and 
p = the correlation coefficient calculated from equation (A-15), 

A closed form solution for this integral is not available, so a MATLAB program is used to 
numerically calculate the probability density flmction. The MATLAB program Puldist.m, 
in Appendix E, calculates the probability density function for the pulse centroid TOA . 

2. Effect of Signal to Noise Ratio on Probability Density 

The plots in Figure A-2 demonstrate the effect of the peak signal to noise ratio on 
the probability density flmction of the pulse TOA. Samples of white Gaussian noise are 
added to the 1.0 microsecond pulse shown in Figure A-1. The variances of the noise are 
chosen to give peak signal to noise ratios of IS, 20, 25, and 30 dB. The variance of the 
noise is calculated from the peak signal to noise ratio with equation (A-22). 




(A-22) 


Where 


V = the variance of the noise required to give the specified peak 
signal to noise ratio for the unit pulse. 
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Figure A-2. Pulse TOA probability density for peak 
signal to noise ratios of 15, 20, 25, and 30 dB 


As expected, as the signal to noise ratio increases, the variance of the error present in the 
estimate of the TOA decreases and the estimate of the pulse TOA becomes more accurate. 

3. Calculated Mean and Standard Deviation of Sampled Pubes 

The mean and standard deviation are calculated for the one microsecond pulse 
shown in Figure A-1. A variety of peak signal to noise ratios and sampling rates are used. 
For all of the sampled pulses the mean time of arrival for the centroid is 0.5 microseconds. 
The calculated standard deviations are shown in Table A-I. 
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Sampling Rate 


S/N 

4MHz 

8 MHz 

12MHz 

l6MHz 

15 dB 

005873 

0 02910 

0 02241 

0 01883 

20 dB 

0 03221 

0 01629 

001258 

001058 

25 dB 

001819 

0 00925 

000713 

000598 

30 dB 

0 01062 

0 00542 

0 00412 

000344 


TABLE A-l. STANDARD DEVIATION FOR 1 0 MICROSECOND 
SAMPLED UInTIT PULSE (microseconds) 


The pulse TOA standard deviations listed m Table A-l indicate that there is a 
tradeoff between sampling rate and peak signal to noise ratio If a lower sampling rate is 
used, a larger peak signal to noise ratio is required to obtain the same TOA error level 
There exists an mversely proportional relationship between the variance of the pulse TOA 
error and the peak signal to noise ratio From Table A>l, the relationship between the 
variance of the TOA error and the peak signal to noise ratio for a sampling rate of 8 MHz 

Van^ = (0.02 9 10)^ ^ ^ _ 

Var2MB (0.00925)2 ' ^ 


B. BURST TOA PROBABILITY DENSITY 

1. Derivation of Probability Density 

The algorithm that calculates the probability density function for the sampled 
pulse is used to calculate the probability density function for a burst of pulses. The 
MATLAB program Burdist.m, listed in Appendix E, generates the burst shown in 
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Figure A-3 and calculates the p ^ility density function for a variety of peak signal to 
noise ratios, PRFs and pulse wiuuis 



Figure A-3. A burst of sampled pulses from a scanning emitter 

The envelope of the burst is approximated as the upper half of a sinewave, and 
the peak signal to noise ratios are calculated using equation (A-22). 

2. Effect of Signal to Noise Ratio on Burst TOA Probability Density 

Figure A-4 demonstrates the effect of signal to noise ratio on probability density 
of the burst TOA. The probability densities for the burst in Figure A-3 are calculated for 
peak signal to noise ratios of 15, 20, 2S, and 30 dB. The noise variances required are 
calculated using equation (A-22). 
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Figure A>4. Burst TOA probability densities for a 2 millsecond burst 
with a pulsewidth of 1.0 miCToseconds, and a PRF of 6 kHz. 

As expected, as the signal to noise ratio is increased, the variance of the error 
present in the estimate of the TOA decreases and the estimate of the burst TOA becomes 
more accurate. 

3. Effect of PRF on the Burst TOA Probability Density 

The probability density is calculated for a burst with a variety of pulse repetition 
frequencies (PRF). The PRFs chosen are 3,6, 9, and 12 kHz. Plots showing the bursts 
used are presented in Figure A-S. 
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PRF = 3kHz 



PRF = 6kHz 




PRF = 9kHz 


PRF= 12 kHz 


Figure A-5. Bursts used to examine the efifect of PRF 
on the TOA probability density 


The other burst characteristics are: 

Burst Length 
Pulsewidth 

Peak Signal to Noise Ratio 
Sampling Rate 


2 0 milliseconds 
1.00 microseconds, 
25 dB, 

8 MHz 
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The calculated probability density functions for the bursts in Figure A-S are plotted in 
Figure A-6. 



Figure A-6. Burst TOA probabibty densities for 
PRFs of 3,6, 9, and 12 kHz 


As the PRF increases the standard deviation of the burst TOA error decreases. 
With an increase in the PRF more energy is contained in the burst and the effect of the 
noise present diminishes 

4. Calculated Mean and Standard Deviation of Sampled Bursts 

The means and standard deviations are calculated for the bursts in Figure A-S for 
a variety of peak signal to noise ratios. The mean time of arrival for the centroid of the 
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burst is 1.0 milliseconds for all of the bursts. The standard deviations are Usted in 
Table A-2. 


Sampling Rate 


S/N 

3000 Hz 

6000 Hz 

9000 Hz 

12000 Hz 

15dB 

0.02392 

0.01665 

0.01353 

0.01185 

20 dB 

0.01337 

0.00925 

0.00755 

0.00655 

25 dB 

0.00754 

0.00524 

0.00428 

0.00373 

30 dB 

0.00428 

0.00300 

0.00248 

0.00218 


TABLE A-2. STANDARD DEVIATION FOR 2.0 MILLISECOND 
SAMPLED UNIT BURST (milliseconds) 


The data in Table A-2 show that the PRFs and the variances of the burst TOA are 
inversely proportional. The relationship between the PRF and the variance of the burst 
TOA is shown below for a signal to noise ratio of 20 dB: 

Variaooo ^ (0 00655)^ - 0 240 « 3000 Hz 
Varjooo (0.01337)2" 12000 Hz 

Similar to the relationship between the PRF and the variance of the burst TOA, the peak 
signal to noise ratio and the variance of the burst TOA are inversely proportional. For a 
PRF of 6000 Hz 


VarndB 

VarjsdB 


(0 01665)2 
(000524)2 


= 10 0s(25dB- I5dB) 
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5. Effect of the Pulsewidth on the Probability Density 


The probability density function of the burst TOA is calculated to er.amine the 
effect of th; pulsewidth on the probability density. A burst is generated with the following 
characteristics; 

Burst Length; 2.0 milliseconds 

Peak Signal to Noise Ratio 20 dB, 

Pulse Repetition Frequency 6000 Hz, 

Sampling Rate 8 MHz 


For this burst the pulse width is varied from 1.0 microseconds to 4.0 
microseconds and the probability density function is calculated and the mean and standard 
deviations are examined. The probability densities for each of the bursts are plotted in 
Figure A-7. 
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Centroid Time of Arrival (milliseconds) 


Figure A-7. Burst TOA probability densities for 
pulse widths of 1.0, 2.0, 3.0, and 4.0 microseconds 


As shown in Figure A-7, the effect of the pulsewidth on the probability density 
function is similar to the effect of the PRF and the peak signal to noise ratio on the 
probability density function. As the pulsewidth increases, the standard deviation of the 
probability density decreases. As the pulse width, and effectively the power present in the 
burst increases, the estimate of the TOA of the centroid of the burst becomes more 
accurate. The pulsewidth and the variance of the TOA of the burst centroid are inversely 
proportional. The relationship between the pulse width and the variance of the probability 
density function is shown below for a signal to noise ratio of 20 dB and a PRF of 6000 
Hz: 
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Var4.o _ (0.00472)^ _ ^ ^ 1.0 ^isec 

Varj.o (0.00926)^ 4.0 ^sec 
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APPENDIX B. LOCI OF CONSTANT TDOA 


A. BURST TDOA PROBLEM 

The burst TDOA for two widely separated receivers being scanned by a constantly 
rotating emitter is directly proportional to the angle formed by the two receivers and the 
emitter. 


y Emitter 



Figure B-1. Angular relationships, burst TDOA problem 
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Using the geometry in Figure B-1, a relationship exists between the bearing and range 


from receiver A to receiver B, the bearing and range to the emitter, and the angle 6. 

_Ra_ Rab 

sin(7i - (0 + 4) - v|/)) sin(6) ^) 

Where SR = the scan rate of the emitter in radians/second, 

Ra = the range to the emitter, 

4* = the bearing to the emitter, 

Rb = the range to receiver B, 

\|/ = the bearing to receiver B, and 
6 = the angle formed by the receivers and the emitter. 

When this equation is rearranged; 


tan(0) = 


sin(4)-\|/) 

SR[^"Cos(4>-m/)] 


(B-2) 


If the assumption is made that the angle 6 is small, then the tangent of 6 is 
approximately equal to the angle 0 and the TDOA can be found from the following 
equation: 


TDO,.,- e _ sin((|i-v) 

SR SR[ft-cos(*-v)] 


(B-3) 


Equation (B-3) is used to plot all of the possible locations that an emitter could be located 
given a TDOA observation. 

B. PULSE TDOA PROBLEM 

A simple relationship exists that can be used to plot the loci of constant TDOA for the 
pulse TDOA problem. The pulse TDOA problem geometry is shown in Figure B-2. 
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Emitter 


y 



Figure B-2. Pulse TDOA problem geometry 

From the law of cosines; 

Rb^ = Ra^+Rab^-2RaRabcos(<|)-\|/) (B-4) 

The difiTerence between Ra and Rb is the distance that a pulse must travel to reach the 
more distant receiver A. The amount of time required to travel this distance is the pulse 
TDOA. This distance D can be found from the following equation; 

D = (TDOA)c (B-5) 

Where TDOA = the pulse time difference of arrival, 

c = the speed of light. 
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Therefore, Rb is: 


Rb = Ra-D 


(B-6) 


Substituting into equation (B>4): 

(Ra - D)2 = Ra^ + Rab^ - 2RaRab cos(4 - v) (B-7) 


(t) = arccos 


' Rab^-2RaD + D^ 
2RaRab 




(B-8) 


Equation (B-8) can be used to plot the possible locations of an emitter given the locations 
of the sensors and the TDOA. 




APPENDIX C. TDOA ORTHOGONALITY 


A. BURST TDOA ORTHOGONALITY 

Equation (C-1) is an approximate reiationsiup used to calculate the loci of constant 
TDOA for the burst TDOA problem 


TDOA = 


sin(^-V) 

SR[^-cos(4>-\(// 


(C-1) 


Where SR = the scan rate of the emitter in radians/second, 

Ra = the range to the emitter, 

4) = the bearing to the emitter, 

Rb = the range to receiver B, 

\|/ = the bearing to receiver B, and 
6 s the angle formed by the receivers and the emitter. 

A measure of the orthogonality of the TDOA observations can be found by taking the 


dot product of the unit vectors tangent to the loci of constant TDOA at the coordinates of 


the emitter. These vectors can be calculated from the slope of the line tangent to the loci 


at the point of interest. The dot product of the unit vectors tangent to the loci will yield 
the cosine of the angle formed by the two loci. The equation for the dot product is given 
below: 


Rt 1 • Rt 2 = cos 0 (C-2) 

Where Rt, = The unit vector tangent to the loci number 1, 

Rt, = the unit vector tangent to loci number 2, and 
0 = the angle formed by the two vectors. 

The cosine of the angle formed by the vectors tangent to the loci is an indication of the 
orthogonallity of the loci. If the cosine of the angle formed by the loci is close to one then 
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the loci are nearl> parallel If the cosine of 6 is zero the loci are orthogonal and lie at right 
an^es to each other The slope, dyt/dxt, of the loci of constant TDOA was found by 
expressing equation (C-1) in Cartesian coordinates and implicitly diflferentiating An 
assumption is made to simplify the problem. If the distance from receiver A to the emitter 
is much larger than the distance to receiver B, the ratio Ra/Rab will be much larger than 
one, and equation (C-I) can be approximated with equation (C-3). 


TDOA = 


Rabsin((j)-v|/) 

SR(Ra) 


(C-3) 


For a given loci the TDOA, Rab, and SR will be constant and are not a function of the 
location of the emitter. The equation can be rewritten in the following form; 


SR(TDOA) sinil)Cos\|/-sinM/cos4) 

Rab ~ Ra ^ 

The following substitutions are made to solve for the equation in terms of the location of 

the emitter xt and yt. 

sin 4) = ^ cos<t> = ^ Ra = [(xt-xa)^+(yt-ya)^]'^ (C-5) 

Ke Rfl 

The resulting equation is; 


SR(TDOA) _ (yt - ya)cos xa)sin \\i 

Rab [(xt - xa)^ + (yt - ya)^] 

The cosv|/ and the siny terms are retained to make the derivations more clear. Implicit 
differentiation was used to calculate the derivative dyt/dxt of equation (C-6). The term on 
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the left side of equation (C-6) is a constant and therefore its derivative is zero. The final 
derivative is given by equation (C-7). 



dxt [(yt - ya)^cos \\i - 2(xt - xaKyt - ya)sin vj/ - (xt - xa)^cos v|/] 


Where xt, yt = the x and y coordinates of the emitter, 

xa, ya = the X and y coordinates of receiver A, and 
V = the bearing fi-om recover A to receiver B. 

The bearing fi’om receiver A to receiver B is calculated fi'om the following equation; 


\^ = arctan 



(C-8) 


The unit vector tangent to the loci of constant TDOA at the coordinates of the emitter is 
given by the following equation; 



Where Rt - the unit vector tangent to the loci, and 

m - the slope of the loci dyt/dxt at the emitter. 


(C-9) 


B. PULSE TDOA ORTHOGONALITY 

The dot product of the unit vectors tangent to the loci of constant TDOA give the 
cosine of the angle between the loci of constant TDOA and a good indication of the 
orthogonality of TDOA observations. The unit vectors tangent to the loci are calculated 
from the slope of the loci. The slope of the loci of constant TDOA cai^ be found by taking 
the derivative dyt/dxt of equation (C-10). 
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";TX)A = 


(xt - xa)^ + (yt - ya)^ 1 - [(xt - xb)^ + (yt - yb)^ 


(C-10) 


Where xt, yt = the x and y coordinates of the emitter, 

xa, ya = the x and y coordinates of sensor A, 

xb, yb = the X and y coordinates of sensor B, and 
c = the speed of li^t. 

An assumption will simplify the problem and yield an approximate, but simpler solution. 

If the distance from the sensors to the emitter is much larger than the distance between the 
sensors, the approaching pulses can be considered plane waves, and the TDOA for the 


pulses arriving at the sensors can be calculated from the geometry in Figure C-1. 
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The TE>OA calculated from the geometry in Figure C~ 1. 


TDOA = 


Rab cos(^ - v) 


(C-11) 


Where Rab = the distance from sensor A to sensor B, 

= the bearing to the emitter, 

\|/ = the bearing from sensor A to sensor B, and 
c = the speed of Ught. 

Using trigonometric identities and the substitutions in equations (C-12) and (C-13), 
equation (C-11) can be rewritten as equation (C'14). 


cos4> = 


_ (xt-xa) _ 

[(xt - xa)^ + 0^ - ya)^] 


(C-12) 


sin4> = 


(yt-ya) 

[(xt - xa)^ + (yt - ya)^] 


(C-13) 


c(TDOA) (xt - xa)cos m/+( yt - ya)sin vy 
Rab " ((xt-xa)2+(yt-ya)2]"^ 


(C-14) 


Implicit differentiation is used to calculate the derivative of equation (C-14). The 
derivative of the left hand side is zero because for a givoi loci the terms on the left hand 
side of equation (C-14) are constant. The equation for the slope of the loci of constant 
TDOA is given in equation (C-15). 


dyt _ (yt - ya)^cos \|/ - (xt - xa)(yt - ya)sin \|/ 
dxt (xt-xa)(yt-ya)cos>|/-(xt-xa)^sinM/ 

Where xt, yt = the x and y coordinates of the emitter, 

xa, ya = the x and y coordinates of sensor A 

xb, yb = the X and y coordinates of sensor B, and 
\|/ = the bearing from sensor A to sensor B. 


(C-15) 
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The bearing from sensor A to sensor B is calculated from the following: 


«|/ = arctan 


[ (yb-ya) ' 

[(xb-xa). 


(C-16) 


The unit vector tangent to the loci of constant TDOA at the coordinates of the emitter is 
given by the following equation: 



(C-17) 


Where 


Rt = the unit vector tangent to the loci, and 
m = the slope of the loci at the emitter. 


APPENDIX D. BURST H** DERIVATION 


The equation that calculates the TDOA observations as a function of the receiver 
locations, the emitter locations, and the scan rate of the emitter, is given by; 


TBOA = ^ 


_ (xt - xa)(yt - yb) - (yt - yaXxt - xb) _ 

[(xt - xa)^ + (yt - ya)^l '^[(xt - xb)^ + (yt - yb)^] 


(D-1) 


Where xt,yt = the x and y coordinate of the emitter, 

xa, ya = the X and y coordinate of receiver A, 

xb, yb = the X and y coordinate of receiver B, and 
SR = the scan rate of the emitter in rad/sec. 

The following MATLAB function was used to evaluate the partial derivative of the 


observation equation (D-1). The terms and derivatives in the MATLAB function were 


calculated using a symbolic processor contained in the software package MATHCAD. 


function [Hx,Hy] = hk3(xt,yt,xa,ya,xb,yb,SR); 

% A function to calculate the derivative of the observation equation 
% This ftmction assumes that all of the receivers are allowed to move, and 
% Their locations are known. The scan rate is known a priori. 

% 

% Richard W. Williamson 
% Dale: 18 February 1994 
% Date Revised: 18 July 1994 
% 

% This program evaluates the partial derivatives of equation (D-1) 


% Define extra variables for clarity. 

xa2 = xa''2; 
ya2 = ya''2; 
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xb2 ® xb''2; 
yb2 * yb''2; 


xt2 = xt^2; 
xt3 = xt'^S; 

yt2 = yt''2; 
yt3 = yt'^S; 

% Denominator of the original TDOA fimction; Equation G^-1) 
G = ((xt-xa)^2 + (yt-ya)'^2 )*( (xt-xb)''2 + (yt-yb)''2 ); 


dgdx = -4*xt*yt*yb + 8*xt*xa*xb - 2*xa*xb2 - 2*xa*yt2 - 2*xa*yb2 
2*xa2*xb - 2*yt2*xb - 4*yt*ya*xt - 2*ya2*xb + 4*xt3 -... 
6*xt2*xb + 2*xt*xb2 + 4*yt2*xt + 2*xt*yb2 - 6*xt2*xa +... 
2*xa2*xt + 2*ya2*xt + 4*xa*yt*yb + 4*yt*ya*xb; 


dgdy = -2*xt2*yb - 4*xt*xa*yt - 2*xa2*yb - 4*yt*xt*xb - 2*ya*xt2 -... 
2*ya*xb2 +8*yt*ya*yb - 2*ya*yb2 - 2*ya2*yb + 4*yt3 +... 

4*xt2*yt + 2*xa2*yt + 2*yt*xb2 - 6*yt2*yb + 2*yt*yb2 -... 

6*yt2*ya + 2*ya2*yt + 4*xt*xa*yb + 4*ya*xt*xb; 

% Numerator of the original TDOA function; Equation (D-1) expanded out 
f = xt*(ya-yb) + yt*(xb-xa) + (xa*yb-xb*ya); 

dfdx = (ya-yb); 
dfdy = (xb-xa); 

% The calculated numerators of the derivative functions 

numdhdx = G*dfdx - 0.5*f*dgdx; 
numdhdy = G*dfdy - 0.5*f*dgdy; 


% Evaluate the partial derivatives 

Hx = (l/SR)*numdhdx/(sqrt(G^3)); 
Hy = (l/SR)*numdhdy/(sqrt(G^3)); 


108 


APPENDIX E. MATLAB PROGRAMS 


% Burst.m 
% 

% Name: Richard W. Williamson 
% Date: 18 Feburary 1994 
% Date Revised: 18 July 1994 
% 

% This program implements an extended Kalman filter to estimate 
% the location of an emitter from the TDOAs of a burst between three 
% receivers. The position of the receivers are known exactly. 

% The algorithm generates noise free TDOA obsenrations and adds white 
% Gaussian noise to the observations. The variance of the white Gaussian 
% noise is specified in the program and can be different for each receiver. 

% Twenty-five observations are processed. 

% 

% The sensors are initially located at the following coordinates: 

% 

% Receiver A (500,100) 

% Receiver B (20000, -2000) 

% Receiver C (-20000,-20000) 

% 

% The TDOA measurements are corrupted by white Gaussian noise. 

% The variance of the noise for all TDOA observations is 1.694e-15. 

% 

% The emitter was located at a range of 500 kilometers from the origin 
% and at a bearing of 90 degrees measured from the x-axis. 

% The coordinates of the emitter were: 

% 

% Emitter: (0, 500) kilometers. 

% 

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 
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% 

% The program outputs four figures. 

% 

% figure(1) The plot of the estimates of the x and y coordinates of the 
% emitter 

% figure(2) The Kalman gains calculated for each observation 
% figure(3) An X-Y plot of the estimates of the coordinates of the emitter. 

% This plot also includes the loci of constant TDOA and the 

% 3a error ellipsoids for the 1st, 8th, 15th, and 22nd estimate 

% of the location. 

% figure(4) A closeup of the steady state estimate of the emitter location. 

% The plot is centered on the steady state estimate of the 

% emitter location and the 3a error ellipsoid corresponding 

% to the steady state estimate is plotted. 

% 

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% The program calls four subroutines. 

% 

% tdoa3.m 
% 

% 

% 

% 

% 

% hk3.m Calculates the derivative of the observation equation given 

the locationof the two receivers, the scan rate, and the 
% estimated location of the emitter. 

% 

% bloci.m Calculates the loci of constant TDOA given the TDOA 

% observation and the locations of the receivers. 

% 

% ellip.m Calculates the error ellipsoids given the location of the point 

% of interest, the error covariance matrix, and the number of 

% standard deviations required. 

% 

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

% 


Calculates the time difference of arrival of a burst 
between two receivers. The locations of the receivers 
are arbitrary. The TDOA is calculated from the locations 
of the receivers, the location of the emitter, and the scan 
rate. 
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clear; 


% Range and bearing to the target 
Rt = 500000; 

Bmgt = 90; 

xt s Rt*cos(Bmgt*pi/180); 
yt =s Rt*sin(Bmgt*pi/180); 

% The number of observations processed. 

k= 1:25; 
len = length(k); 

% The coordinates of receivers A, B, and C as a function of k 

xa = 500*ones(1,len) + 00*k;; 

ya = 100*ones(1,len) + 30*k;; 

xb = 20000*ones(1.len) + 30*k; 

yb = -2000*ones(1.len) - 30*k; 

xc = -20000*ones(1,len) - 30*k; 

yc = "20000*ones(1 ,len) - 30*k; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 

% Initializing variables 

Po= 1.00e+20*eye(2); 

Q = 4.0e+06*eye(2); 

RKab= 5.660e-09; 

RKac= 5.660e-09; 

RKbc= 5.660e-09; 

Srate = 2*pi; 

X = zeros(2,len); 

Pk = zeros{2,2); 

HK = zeros(3,2); 
g = zeros(2,len); 
z_tru = zeros(3,len); 

Z = zeros(3.len); 

Zhat = zeros(3,len); 


% The initial error covariance matrix 
% Variance of the state excitation noise. 

% The variance on the TDOA error Rcvr A-B 
% The variance on the TDOA error Rcvr A-C 
% The variance on the TDOA error Rcvr A-C 
% The scan rate of the emitter 

% The state i iatrix 
% The error covariance matrix 
% The linearized derivative of h(k) 

% The Kalman gains 
% The vector of true observations 
% The vector of noisy observations 
% The vector of estimated obsen/ations 
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% Calculation of the observations noise free observations 
fork* 1:len, 

z_tru(1,k) = tdoa3(xt,yt,xa(k),ya(k),xb(k),yb(k),Srate); 
2_tru(2,k) = tdoa3(xt,yt,xa(k).ya(k),xc{k),yc(k),Srate); 
z_tru(3,k) = tdoa3(xt,yt,xb(k),yb(k),xc(k),yc(k).Srate); 
end; 

% Calculate the noise on the TDOA estimates 

rand(’normar): 

n = sqrt(RKab)*rand(3,len); 

% Form the matrix of noisy observations. 

Z = 2 _tru + n; 

% A priori estimate of the location of the emitter 
Xo = [0 10000]': 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Initialize the Extended Kalman filter 

Pk = Po: % Initialize the Pk the error covariance matrix. 

X(:,1) = Xo; % Initialize the state vector 

fork * 1:len, 

k 

% Processes the observation for receiver A and B 

Calculates the derivative of the Rcvr A-B TDOA equation wrt x and y 
[hlab h2ab] = hk3(X(1.k).X(2.k).xa(k).ya(k).xb(k),yb(k),Srate): 

% Forms the Hk matrix 
HKab = [hlab h2ab]; 

% Calculates the next PK based upon the previous PK 


Pk = Pk + Q: 
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% Calculates the kalman gains 

GKab = Pk*HKab'*jnv( HKab*Pk*HKab’ + RKab ); 

% Calculates the estimate of the observation based upon the 
% estimate of the states. 

Zhat(1,k) = tdoa3(X(1.k).X(2.k),xa(k).ya(k).xb(k).yb(k).Srate): 

% Calculates an update of the state based upon the difference in 
% the actual observation and the estimate of the observation. 
X(:,k) = X(:,k) + GKab*( Z(1.k)-Zhat(1.k)); 

% Updates the variance matrix Pk based upon the observation 
Pk = (eye(2) - GKab*HKab)*Pk; 

% Processes the observation for receiver A and C 

% Calculates the derivative of the TDOA equation wrt x and y 
[h1ac h2ac] hk3(X(1.k).X(2.k).xa(k).ya(k).xc(k).yc(k),Srate): 

% Forms the Hk matrix 
HKac = [hlac h2ac]: 

% Calculates the next PK based upon the previous PK 
Pk = Pk + Q: 

% Calculates the kalman gains 

GKac = Pk*HKac'*inv( HKac*Pk*HKac’ + RKac); 

% Calculates the estimate of the observation based upon the 
% estimate of the states. 

Zhat(2.k) = tdoa3(X(1,k),X(2,k).xa(k).ya(k).xc(k).yc(k).Srate); 

% Calculates an update of the state based upon the difference in 
% the actual observation and the estimate of the observation. 
X(:.k) = X(;,k) + GKac*( Z(2.k)-Zhat(2.k)); 

% Updates the variance matrix Pk based upon the observation 
Pk = (eye(2) - GKac*HKac)*Pk: 
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% Processes the observation for receiver B and C 

% Calculates the derivative of the TDOA equation wrt x and y 
[h1bc h2bc] = hk3(X(1,k),X(2,k),xb(k).yb(k).xc(k).yc(k),Srate): 

% Forms the Hk matrix 
HKbc = [h1bc h2bc]: 

% Calculates the next PK based upon the previous PK 
Pk = Pk + Q: 


% Calculates the kalman gains 

GKbc = Pk*HKbc“inv( HKbc*Pk*HKbc* + RKbc); 

% Calculates the estimate of the observation based upon the 
% estimate of the states. 

Zhat(3,k) = tdoa3(X(1,k),X(2,k).xb(k),yb(k),xc(k),yc(k),Srate): 

% Calculates an update of the state based upon the difference in 
% the actual observation and the estimate of the observation. 
X(:.k) = X(:.k) + GKbc*( Z(3,k)-Zhat(3.k)); 

% Updates the variance matrix Pk based upon the observation 
Pk = (eye(2) - GKbc*HKbc)*Pk: 

% Projects the next state based upon the current state 

X(;.k+1) = X(:.k): 

% Save the error covariance matricies 
ifk== 1 
P = Pk; 
else, 

P = IPiPk]: 
end; 

% Save the Kalman gains 
g(;.k) = GKbc; 


end; 


% End of the filtering routine 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%o^%o^% 
% Output the filtered estimates of x and y 


t = 1 :len; 
figure(1) 

plot(t.X(1 .t)/1000.t.X(2.t)/1000) 
axis([1 25 -50 600]) 
pause; 

% Plot of the Kalman gains 

figure(2) 

cig 

plot(t.g(1.:).t.g(2.:)) 
title('Kalman Gains') 
xlabeK'Number of Observations ) 
ylabel('Magnitude') 
pause 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Calculate the steady state estimate of emitter location 

xtavg = mean(X(1,k-10:k)); 
ytavg s mean(X(2.k-10:k)); 


% Calculate the coordinates of the loci of constant TDOA based 
% upon noise free TDOA observations 

Tab = bloci(xa(k),ya(k).xb(k).yb(k).z_tru(1.k),Srate); 

Tac = bloci(xa(k),ya(k).xc(k),yc(k),z_tru(2.k),Srate); 

Tbc = bloci(xb(k),yb(k),xc(k),yc(k).z_tru(3,k),Srate); 


% Plot the locus of emitter locations and actual emitter locations 

figure(3) 

cig 

% Plots coordinates of the estimates of the emitter location 
% axis([-500 500 -30 1000]); 
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plot(X(1 ,:)/1000,X(2,:)/1 OOO.r) 
hold on 


% Plots the actual emitter location 
plot(xt/1000.yt/1000.'*g') 

% Plots the final locations of the receivers 
plot(xa(k)/1000.ya(k)/1000.”*g*) 
plot(xb(k)/1000,yb(k)/1000,‘‘g’) 
plot(xc(k)/1000.yc(k)/1000.’*g') 

% Plots the locus of possible emitter locations 
plot( Tab(1.:)/1000, Tab(2,:)/1000.’:b*) 
plot( Tac(1,:)/1000. Tac(2.:)/1000.':b’) 
plot( Tbc(1.:)/1000. Tbc(2.:)/1000.':b') 

% Plot three sigma error elipsoids for various points 

C = 3.00; % The number of standard deviations plotted in error ellipsoids 

fork= 1:7:len"1, 

P1 * [P(2*k-1,1) P(2*k-1,2); P(2*k,1) P(2*k,2)l: % The Pk Matrix 

lxg,y1 ,y2,Emax,Emin] = ellip(P1 ,C,X(1 ,k),X(2,k)): % Calculate ellipse 

plot((xg)/1000,(y1 )/1000,'g-’,(xg)/1000.(y2)/1000,*g-’) % Plot the ellipse 

end; 

% axis; 
hold off 

axis([-250 250 -50 600]); 
pause 

% Plot the locus of emitter locations and actual emitter locations 

figure(4) 

cig 

% Plots coordinates of the estimates of the emitter location 

plot(X(1 ,:)/1000,X(2,:)/1000,'r’) 
hold on 
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% Plots the actual emitter location 
plot(xt/1000.yt/1000;*g'): 

% Plots the locus of possible emitter locations 
plot( Tab(1.:)/1000. Tab(2.:)/1000.‘:b'): 
plot( Tac(1.:)/1000, Tac(2.:)/1000.':b'); 
plot( Tbc(1.:)/1000, Tbc(2.:)/1000.':b'): 

% Plots the error ellipsoids 

[xg,y1,y2,Emax.Emin] = ellip(Pk.C.xtavg.ytavg); 

plot((xg)/1000.(y1)/1000;g-'.(xg)/1000.(y2)/1000;g-') 

Major = sprintf('%5.1f,Emax/1000): 

Minor = sprintf('%5.1f,Emin/1000); 

text((xt-0.50*Emax)/1000,(yt+0.85*Emax)/1000,fMajor Axis = ’ Major' km*]) 
text((xt-0.50*Emax)/1000,(yt+0.75*Emax)/1000,[Minor Axis = * Minor ’ km*]) 

% axis; 
hold off 

axis([(xt-Emax)/1000 (xt+Emax)/1000 (yt-Emax)/1000 (yt+Emax)/1000]); 

% axis('square') 

% 

% 

% 
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% tdoa3.m 


% 

% Name; Richard W. Williamson 
% Date: 18 Feburary 1994 
% Date Revised; 18 July 1994 
% 

% This program calculates the burst observation equation h(x(k)). 


% 

% Input; The coordinates of receiver A (meters); xa, ya 

% The coordinates of receiver B (meters); xb, yb 

% lire coordinates of the emitter (meters); xt, ^ 

% Scan rate of the emitter: SR 

% 

% Output The burst time difference of arrival for 
% the specified receiver loactions. T 

% 

% 

function T = TDOA3(xt,yt,xa,ya,xb,yb,SR) 


% Numerator of the observation equation 
f = xt*(ya-yb) + yt*(xb-xa) + (xa*yb-xb*ya); 

% Denominator of the observation equation 
G = ((xt-xa)''2 + (yt-ya)^2 )*( (xt-xb)''2 + (yt-yb)'^2); 

% Calculate the TDOA 
T = f/(SR*sqrt(G)); 

% 


118 




% HK3.m 
% 

% Name: Richard W. Williamson 
% Date: 18 Feburary 1994 
% Date Revised: 18 July 1994 
% 


% This program calculates the derivative of the observation equation h(x(k)). 
% 


% Input: 
% 

% 

% 

% 

% Output 
% 

% 


The coordinates of receiver A (meters): xa, ya 

The coordinates of receiver B (meters): xb, yb 

The coordinates of the emitter (meters): xt, ^ 

Scan rate of the emitter: SR 

The vector of partial derivaties of h(x(k)) 

with respect to xt and yt. [Hx.Hy] 


function [Hx.Hy] = hk3(xt,yt,xa.ya,xb,yb,SR); 


% A function to calculate the derivative of the 
% H matrix. This function does not assume that one of 
% The receivers is fixed with respect to the stationary 
% coordinate system. This assumes that the scan rate is known. 
% 

% This program is an implementation of equation ( ) 

% 

% Define variable for easier understanding 


xa2 = xa''2; 
ya2 = ya''2; 


xb2 = xb''2: 
yb2 = yb'‘2: 

xt2 = xt''2: 
xt3 = xt''3; 

yt2 = yt''2; 
yt3 = yt''3; 
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% 

% The original TDOA function: 


G = ((xt-xa)''2 + (yt-ya)''2 )*( (xt-xb)^2 + (yt-yb)''2); 

dgdx * -4*xt*yt*yb + 8*xt*xa*xb - 2*xa*xb2 - 2*xa*yt2 - 2*xa*yb2 
2*xa2*xb - 2*yt2*xb - 4*yt*ya*xt - 2*ya2*xb + 4*xt3 -... 
6*xt2*xb + 2*xt*xb2 + 4*yt2*xt ^ 2*xt*yb2 - 6*xt2*xa +... 
2*xa2*xt + 2*ya2*xt + 4*xa*yt*yb + 4*yt*ya*xb: 

dgdy = -2*xt2*yb - 4*xt*xa*yt - 2*xa2*yb - 4*yt*xt*xb - 2*ya*xt2 - 
2*ya*xb2 +8*yt*ya*yb - 2*ya*yb2 - 2*ya2*yb + 4*yt3 +... 
4*xt2*yt + 2*xa2*yt + 2*yt*xb2 - 6*yt2*yb + 2*yt*yb2 -... 
6*yt2*ya + 2*ya2*yt + 4*xt*xa*yb + 4*ya*xt*xb; 

f = xt*(ya-yb) + yt*(xb-xa) + (xa*yb-xb*ya); 


dfdx = (ya-yb): 
dfdy = (xb-xa); 

% The calculated numerators of the derivative functions 

numdhdx = G*dfdx - O.STdgdx; 
numdhdy = G*dfdy - O.STdgdy; 

% The actual derivative functions 

Hx = (1/SR)*numdhdx/(sqrt(G''3)); 

Hy = (1/SR)*numdhdy/(sqrt(G''3)): 
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Bloci.m 


% 


% Name: Richard W. Williamson 
% Date: 18 Feburary 1994 
% Date Revised: 18 July 1994 
% 


% 

% 

% 

% 


% 

% 

% 

% 

% 

% 


This program calculates the coordinates of the loci of constant TDOA for 
a pair of receivers. 


Input: The coordinates of receiver A (meters): xa, ya 

The coordinates of receiver B (meters): xb, yb 

The Time Difference of Arrival (sec): TD 

Scan rate of the emitter: SR 

Ou^ut The vector of x and y coordinates of the loci 

of constant TDOA. T 


% 

function T = bioci(xa,va,xb,yb,TD,Srate) 


% Polar plot of the receivers, the emitter,and the locus of possible 
% emitter location. 

% The angles of phi to be plotted, 
phi > 0:pi/100:pi; 


% Calculates the range and bearing to the receivers. 
Rab = sqrt( (xb-xa)''2 + (yb-ya)''2 ); 

psi = atan2(yb-ya,xb-xa)*ones(1,1:length(phi)); 

% initialize the range variables 
Ra = zeros(1,length(phi)), 

Ra = Rab*( sin(phi-psi)/tan(TD*Srate) + cos(phi-psi)); 

dxa = xa*ones(1,1:length(phi)); 
dya = ya*ones(1,1:length(phi)): 

% Plots the locus of possible emitter location. 

T = [Ra.*cos(phi)+dxa ; Ra.*sin(phi)+dya J 
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% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 


pulse.m 

This program implements an extended Kalman filter to estimate 
the location of an emitter from the TDOA of a pulse between two 
sensors. Two receiver platforms are used and the their location is 
known exactly. Each receiver platform consists of two sensors The 
algorithm generates noise free TDOA observations for the known receiver 
and emitter locations and adds white Gaussian noise to the observations. 
The variance of the white Gaussian noise is specified in the program 
and can be different for each receiver platform. Sixty obsevations 
are processed. 


The sensors are located at the following coordinates: 


Receiver 1 
Receiver 1 
Receiver 2 
Receiver 2 


Sensor A (0,0) 
Sensor B (500,0) 
Sensor A (20000,0) 
Sensor B (20500,0) 


The TDOA measurements are corrupted by white Gaussian noise. 
The variance of the noise for both receivers is 1.694e"15. 


The emitter was located at a range of 30 kilometers from the origin 
and at a bearing of 60 degrees measured from the x-axis. 

The coordinates of the emitter were: 


Emitter: (15, 26) kilometers. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
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% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 


The program outputs four figures. 

figure(1) The plot of the estimates of the x and y coordinates of 
the emitter 

figure(2) The Kalman gains calculated for each observation 

figure(3) An X-Y plot of the estimates of the coordinates of 

the emitter. 

This plot also includes the loci of constant TDOA and the 
3a error ellipsoids for the first, 20th, 39th, and 58th 
estimate of the emitter. 

figure(4) A closeup of the steady state estimate of the emitter 
location. The plot is centered on the steady state 
estimate of the emitter location and the 3a error 
ellipsoid corresponding to the steady state estimate 
is plotted. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


The program calls two subroutines. 

ploci.m Calculates the loci of constant TDOA given the TDOA 
observation and the locations of the receivers, 
ellip.m Calculates the error ellipsoids given the location of the point 
of interest, the error covariance matrix, and the number of 
standard deviations required. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


Initialize the variables 


c = 3.00e+08; 

R1 = 1.694e-15: 

R2 = 1.694e-15; 

Q = 5.0e+02*eye(2): 
PO = 5.0e+20*eye(2): 


% The speed of light 

% The covariance of measurement error receiver 1 
% The covariance of measurement error receiver 2 
% Variance of Plant excitation noise 
% Initial error covariance 


% Choose to calculate the error covariance of 3 sigma 
C = 3.0, 


t = (0:60); % Time vector total of 61 observations 

I = ones(1 ,len); % occurring 50 micro seconds apart 
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% Initialize storage for Kalman gains and the states, 
g = zeros(2,length(t)); 

X = zeros(2,len): 

% Emitter location 

Rt - 30000; % Range to the emitter 

Bt = 60.00; % Bearing to the emitter from the origin 

xt = Rt*cos(Bt*pi/180); % X coordinate of the emitter 

yt = Rt*sin(Bt*pi/180); % Y coordinte of the emitter 

% A priori estimate of the location of the target; 

XO = (10000;50001; 

% 

% The location of the emitter specified for all point in time 

xa1 = 0*l + 0*t; % Receiver 1 Sensor A moves in a straight line 

ya1 = 0*l + 0*t; % at a constant speed 

xb1 = 500*l + 0*t; % Receiver 1 Sensor B moves in a straight line 

yb1 = 0*l + 0*t; % at a constant speed 

xa2 = 20000*1 + 0*t; % Receiver 2 Sensor A moves in a straight line 

ya2 = 000.00*1 + 0*t; % at a constant speed 

xb2 = 20500*1 + 0*t; % Receiver 2 Sensor B moves in a straight line 

yb2 = 000*1 + 0*t; % at a constant speed 


% Calculate the TDOA observations for the two receivers 
% Equation 

% z(k+llk) 

% Calculate the observations for receiver 1 

Z1 = (1/c)*( sqrt( (xt-xa1).''2 + (yt-ya1).''2 ) -... 
sqrt( (xt-xb1).''2 + (>1-yb1).''2 )); 


% Calculate the observations for receiver 2 
Z2 = (t/c)*( sqrt( (xt-xa2).''2 + (yt-ya2).''2 ) -... 


(^) - xa(k))^ + (>t(k) - ya(k))^l'^ - F (xt(k) - xb(k))^ + (yt(k) - yb(k))^ 
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sqrt( (xt-xb2).''2 + (yt-yb2).^2 )); 


% Inject zero mean, noise into the observations 

rand('normar): 

n1 = sqrt(R1)*rand(1,len); 

n2 = sqrt(R2)*rand(1,len); 

Z1n = Z1 +n1: 
iZ2n = Z2 + n2: 

% initialize the Kalman filter 

X(;, 1) := XO; % A priori estimate of the emitter location 

PK = PO; % Initial error covariance matrix 

% Calculate the filtered estimates of the emitter location 

fork = 1:len-1. 

% Process TDOA observation Receiver #1 

/O /O A)/O /OyO/O/il/OAl/O/DAt/O/O/O/O /«/O /D /O Al Ai /O A> 70^ A>A>AlA>A>A>A>A>A) 


% Calculate the HK matrix for receiver platform 1 

% Equation H/Tf'i -I U « (^)-xfa) "j if _ (^)-yb) "\ 

Ra(k) Rb(k) J 'V Ra(k) Rb(k) J 

% Equation Ra(k) = [(xt(k)-xa)^+ (yt(k)-ya)^]^^ 

% Equation Rb(k) = [(Jrt(k) xb)^ + (yt(k) - yby ] 

Ra1 = sqrt( (X(1,k)-xa1(k)).''2 + (X(2,k)-ya1(k)).''2 ); 

Rb1 = sqrt( (X(1,k)-xb1{k)).''2 + (X(2.k)-yb1{k)).''2 ); 

h1x = (1/cr( (X(1,k)-xa1(k))./Ra1 - (X(1,k)-xb1(k))./Rb1 ); 
h1y = (1/c)*( (X(2,k)-ya1(k))./Ra1 - (X(2,k)-yb1(k))./Rb1 ); 


HK = [h1xh1y]; 


% Calculate the estimates of the TDOA based upon the estimate 
% of emitter location and location of receiver platform 1 
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% Equation 
% z(k+llk) 

Zhat1 = (1/c)*(sqrt( (X(1,k)-xa1(k))''2 + (X(2.k)-ya1(k))''2 ) -... 
sqrt( (X(1.k)-xb1(k)r2 + (X(2.k).yb1(k)r2 )); 

% The new error covariance is the same as the old plus Q. 

PK = PK Q; 

% Calculate the Kalman gains 
GK = PK*HK'*inv(HK*PK’‘HK + R1); 

% Calculate the updated error covariance matrix 
PK = (eye(2) - GK*HK)*PK; 

% Calculate a smoothed estimate of the bearing to the emitter 
X(;,k) = X(:,k) + GK*(Z1n(k)-Zhat1): 

% estimate of the emitter location for receiver #2 is same as previous 

X(;,k) = X(;.k): 

% Process TDOA observation receiver #2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% Calculate the HK matrix for receiver platform 2 

„/ ,, .. TT/i-N r iTcxtdchxa) (3h(k)-xb)^ 1 (^(k)-yb)'\ 

% Equation H(k) = -- 

% Equation Ra(k) = [(xt(k) - xa)^ + (^(k) - ya)^ J ^ 

% Equation Rb(k) = [(xt(k) - xb)^ + (^(k) - yb)^ ] ^ 

Ra2 = sqrt( (X(1,k)-xa2(k)).''2 + (X(2.k)-ya2(k)).''2 ); 

Rb2 = sqrt( (X(1,k)-xb2(k)).^2 + (X{2,k)-yb2(k)).''2 ); 
h2x = (1/c)*( (X(1,k)-xa2(k))./Ra2 - (X(1,k)-xb2(k))./Rb2 ); 
h2y = (1/c)*( (X(2,k)-ya2(k))./Ra2 - (X(2.k)-yb2(k))./Rb2 ); 


(^) - xa(k))^ + (^) - ya(k))* 1 - f (xt(k) - xb(k))^ + (^) - yb(k))^ 


126 





HK = [h2x h2y]; 


% Calculate the estimates of the TDOA based upon the estimate 
% of the emitter location and location of receiver Platform 2 


% 


% 


Equation 

z(k + 1 Ik) = 


(^) - xa(k))^ + (yt(k) - ya(k))^ 


- (^) - xb(k))^ + (^) - yb(k))n 
^ - 


Zhat2 = (1/c)*(sqrt( (X(1.k)-xa2(k))''2 + (X{2.k)-ya2(k))'^2 ) -... 
sqrt( (X(1.k)-xb2(k))'^2 + (X(2.k)-yb2(k))>'2 )); 

% The new error covariance is the same as the old plus Q. 
PK = PK + Q; 


% Calculate the Kalman gains 
GK = PK*HK*inv(HK*PK*HK’ + R2); 

% Calculate the updated error covariance matrix 

PK = (eye(2) - GK*HK)*PK: 

% Calculate a smoothed estimate of the bearing to the emitter 
X(:.k) = X(:,k) + GK*(Z2n(k)-Zhat2): 

% Estimate of the location for next observations is same as previous 
X(;,k+1) = X(:,k); 

% Save the error covariance matricies 
if k== 1 
P = PK; 
else, 

P = [P;PK]; 
end; 

% Save the Kalman gains 
g(:.k) = GK; 

end; % end of the Kalman filtering routine 
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% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 

% 

% Plot the output 
% 

% 


% Plot of the estimate of the x and y coordinates 

figure(1) 

cig; 

plot(t,X(1 ,;)/1000.t.X(2.:)/1000) 
titleC Estimate of x and y coordinates'); 
xlabei('Number of Observations') 
ylabel('X location in kilometers') 
pause 

% Plot of the Kalman gains 

figure(2) 

cig 

plot(t.g(1.:).t.g(2.:)) 
title('Kalman Gains') 
xlabeK'Number of Observations') 
yiabel('Magnitude') 
pause 

% Calculate the steady state estimate of emitter location 
xtavg = mean(X(1,k-10:k)); 
ytavg = mean(X(2,k-10:k)); 

% Calculate the loci of constant TDOA based upon noise free TDOA 
% observations 

R1a = 5000:5000:40000; 

% T1 is a vector of the X and Y coordinates of the loci. 

T1 = ploci(xa1(k),ya1(k),xb1(k),yb1(k).Z1(k).R1a); 

R2a = 5000:5000:40000; 

% T2 is a vector of the X and Y coordinates of the loci. 

T2 = ploci(xa2(k).ya2(k),xb2{k).yb2(k),Z2(k).R2a); 
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% X-Y Plot of the estimate of the coordinates of the emitter 
figure(3) 


% True emitter location 

subplot(111). plot(xt/1 OOO.yt/1 OOO,-) 

hold on 

% Plots the final locations of the sensors 

plot(xa1 (k)/1000,ya1 (k)/1000;r*'.xb1 (k)/1000.yb1 (k)/1 OOO/r*') 

plot(xa2(k)/1000.ya2(k)/1000.'r*'.xb2(k)/1000.yb2(k)/1 OOO/r*') 

% A priori estimate of emitter location 
plot(X0(1 )/1000,X0(2)/1 000 ,'+') 

% Estimates of emitter location 
plot(X(1 .:)/1000.X(2.:)/1000) 

% Plot three sigma error elipsoids for various points 
for k = 1:19:len-1, 

% Selects the error covariance from those saved in vector P. 
P1 = [P(2*k-1,1) P(2*k-1,2): P(2*k,1) P(2*k,2)];' 

% Calculates the coord, of ellipse. 

(xg,y1,y2,Emax,Emin] = ellip(P1,C,X(1,k),X{2,k)); 

% Plots on graph 

plot((xg)/1000.(y1 )/1000.'g-',(xg)/1000.(y2)/1000,'g-') 

% Plots the loci of constant TDOA 
plot(T1 (1 ,:)/1000,T1 (2,:)/1000,':') 
plot(T2(1,:)/1000,T2(2,:)/1000,’:’) 

hold off 

title('Estimate of Position of the Emitter: tgt') 

% axis; 

axis([0 30 0 30]); 
pause; 


end 
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% Plot the closeup of the estimate of the emitter location 
figure(4) 

% True Emitter location 

subplot(111). plot(xt/1000,yt/1 OOO,'*') 

hold on 

% Estirr ? of emitter location 
plot(X( 000.X(2,:)/1000); 

% Calculate and plot error ellipsoids for steady state estimate 

[xg,y1,y2,Emax.Emin] = ellip(PK,C.xtavg.ytavg): 
plot((xg)/1000,(y1)/1000.'g-'.(xg)/1000.(y2)/1000,'g-‘) 

% Plot on the graph the major and minor axis. 

Major = sprintf('%5.1f.Emax/1000); 

Minor = sprintf('%5.1f.Emin/1000); 

text((xt-0.50*Emax)/1000,(yt+0.85*Emax)/1000,[’Major Axis =' Major' km']) 
text{(xt-0.50*Emax)/1000,(yt+0.75*Emax)/1000,fMinor Axis =' Minor' km*]) 

% Plot the loci of constant TDOA 
plot(T1 (1 ,:)/1000,T1 (2,:)/1000,':') 
plot(T2(1,:)/1000,T2(2,;)/1000,’:') 

hold off 

titleCCIoseup of Position of the Emitter;’) 

% axis; 

axis{[{xt-Emax)/1000 (xt+Emax)/1000 (yt-Emax)/1000 (yt+Emax)/T000]); 
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% Ploci.m 
% 

% This program calculates the coordinates of the loci of constant TDOA for 


% a pair of sensors. 

% 

% Input: The coordinates of sensor A (meters): 

xa.ya 

% 

The coordinates of sensor B (meters): 

xb, yb 

% 

The Time Difference of Arrival (sec): 

TDOA 

% 

% 

The range from sensor A over which 
the loci of constant TDOA will be plotted: 

R 

% 

% Output 
% 

The vector of x and y coordinates of the loci 
of constant TDOA. 

T 

% 

% 

function T = 

ploci(xa,ya,xb.yb,TDOA,R); 



% Calculate the distance from sensor A to Sensor B 
Rab = ones(1,length(R))*sqrt((xb-xa)''2 + (yb-ya)'^2): 


% Calculate the bearing of sensor B from sensor A 
psi * ones(1,length(R))*atan2(yb-ya,xb-xa); 


% Calculate the difference in path length for the pulse 
D = 3.0e+08*TDOA*ones(1.length(R)); 


Ecjuotion (|> = arccos 


[(R-D)^-R^-R^] 


-2R,bR 


% Calculate the bearing required for each range 

phi = acos( ((R-D)/2 - (R).''2 - Rab.''2 )./(-2*R.*Rab )) + psi; 


J = ones(1 ,length(R)); 


% Calculate the x and y coordinates 
T = [R.*cos(phi)+xa*J: R.*sin(phi)+ya*JJ; 
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% ellip.m 
% 

% This program calculates the coordinates of the 2D error ellipsoids given the 
% error covariance matrix and the location of the estimate. 

% 


% 

input; 

The 2D Error Covariance Matrix 

PK 

% 


The number of standard deviations 

C 

% 

% 

% 


The X and y coordinates of the estimate 

xt. yt 

Output 

The X vector along which the ellipsoid is 


% 

% 


plotted. 

The y vectors corresponding to the upper 

xgout 

% 


and lower parts of the ellipsoid 

y1out, y2out 

% 


The length of the major axis of the ellipse 

Axmax 

% 


The length of the minor axis of the ellipse 

Axmin 


% 

function [xgout,y1out,y2out,Axmax,Axmin] = ellip(PK.C,xt,yt) 


% Inputs the error covariance matrix, the number of 
% standard deviations, and the location where plotted 

% Calculate the inverse of the error covariance matrix 
P = inv(PK): 

P11 = P(1.1); 

P12 = P(1.2): 

P21 = P(2,1); 

P22 = P(2,2): 

% Calculate the upper and lower limits on the x axis 

% Equation 
f=sqrt(P11*C''2): 

% Define the vector along the x axis over which the ellipsoid lies 
xg = -f:(2*f)/1000:f; 
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% Calculate the upper ar^d lower curves of the ellipse. 


% Equation 



(PjlC^-JCtO 


y1 = -(P12/P11)*xg + sqrt( ((PirP22 - P12''2)/Pir2)*... 
(P11*C*2*ones(1 ,length(xg)) - xg.''2) ); 

y2 = -(P12/P11)*xg - sqrt( ((P11*P22 - P12''2)/P11''2)*... 

(P11*C''2*ones(1 ,length(xg)) - xg.'^2) ); 

% Calculate the distance of the ellipse from the center 

Rellip =5 sqrt( xg.''2 + y1 .''2 ); 

% Calculate the major and minor axis of the ellipse 

Axmax = abs(2*max(Rellip)); 

Axmin = abs(2*min(Reliip)); 

% Calculate the coordinates of the ellipse over the estimate 

y1out = y1+yt*ones(1,length(xg)); 

y2out = y2+yt*ones(1,length(xg)); 

xgout = xg+xt*ones(1,length(xg)); 

% 
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% Burdist.m 
% 

% Name; Richard W. Wiliiamson 
% Date: 18 February 1994 
% Date Revised; 18 February 1994 
% 

% This program generates a simulated burst of pulses generated by a 
% circularly scanning emitter. The length of the burst is 2.0 milliseconds the 
% maximum amplitude is 1.0, and the PRF varies from 3 kHz to 12 kHz. 

% The program calculates the probability density function for the centroid of 
% the burst for a peak signal to noise ratio of 20 dB. 

% 

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 

% The program outputs one figure. 

% 

% figure(1) A plot of the probability density function for the burst centroid. 

% 

% o^o^o^o/o^o^o^o^o^ozo^o^ozo/,o^o^o^q/qzo^o^o^o^o^ozo^Q,^o^o^o^ozo^ozo^ozq/q/ 
/0/vAl/DA)/0/OAi/OAI/O/O/D/OAlA)/O/0/D/v/D/D/D/D/D /D/D /0/0/0f0A>/0/O/0f0/0 

% 

% Input Variables; 

% 

% BW = Beam width in degrees 
% SR = Scan rate of emitter in degrees/sec. 

% PRF = Pulse Repetition Frequency in Pulses/sec 
% PW = Pulse width in seconds. 

% FS = Sampling Frequency in samples/sec 
% 

% 

% Output Variables; 

% 

% Bur = The noise free Burst vector 
% 

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
clear; 

% Initialize the variables 

BW= 1.00; 

SR = 500.0; 

PRFk = [3000 6000 9000 12000] 

PW= 1.0e-06; 

FS = 8.0e+06; 
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muz = {000 OJ; 
sigz = [000 0], 


zk = {0.90:(1.10-0.900)/200:1.10)/1000; 
dtz = zk(2)-zk(1); 

for m = 1:length(PRFk) 

PRF = PRFk(m); 

% Calculated Variables 


Bur_Len = BW/SR; % Burst length in seconds; 

N_Bur = Bur_Len*FS; % Number of samples per burst 


N_PW = PWTS; % Number of samples per pulse width 

N_PRF = (1/PRF)*FS: % Number of samples per PRF 


Bur = zeros(1 ,N_Bur); % Initialize the Burst vector 
t = 0:N_Bur-1; % Initialize the time vector 


% Load the individual pulses in the burst vector. 

forl = 1:N PRF;N Bur, 

Bur(l:l+N_PW-1)‘= ones(1,N_PW); 
end; 


% Define burst envelope and form noise free burst vector 

Bur_env = 1.0*sin(pi*l/N_Bur); 

Bur = Bur.*Bur_env; 


dt = 1/FS; 

% The variance of the noise power added to the pulse 
V = 0.010; 

ns = Bur; 

N = length(ns); % Length of the burst. 
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% Numerator = x Denominator = y; 

% Initialization of the mean and sigma variables 

mux = 0; 

muy = 0; 

sigx = 0; 

sigy = 0; 

cov = 0; 


% Calculate the statistics of the numerator and 
% the denominatior of the centroid equation. 


Np = 0: 


fork= 1:N, 

Np = Np + 1: 
if ns(k) > 0; 

mux = mux + k*dt*ns(k); 
sigx = sigx + ((k*dt)''2)*V; 
muy = muy + ns(k); 
sigy = sigy + V; 
cov = cov + k*dt*V; 

end; 

end; 

sigx = sqrt(sigx); 
sigy = sqrt(sigy): 

rho = cov/(sigx*sigy) % The correlation coefficient 

c_bar = mux/muy 


dty = (10*sigy)/60; 

yk = (muy-5*sigy):dty:(muy+5*sigy); 

p = zeros(length(yk),length( 2 k)); 


136 



% Calculate the probability distribution for f(zy.y) 
for I = 1 :length(zk), 

for k = 1 ;length(yk), 

y = yk(k): 
z = zk(l): 

g1 = (y*z-mux)/sigx; 
h1 = (y-muy)/sigy; 

fxy = (1/(2*pi*sigx*sigy*sqrt(1-rho''2)))*... 

exp( (-1/(2*(1-rho''2)) )*(g1''2 -2*rho‘g1*h1 + h1''2)); 

p(k.l) = fxy; 

end; 

end; 


% Initialize the variables, 
z = zeros(1 ,length(zk)); 

% Integrate the probability distribution across y to 
% calculate the marginal density for z the centroid variable. 

for I = 1:length(zk), 

H = 0; 

for k = 1:length(yk)-1, 

H = H + dty*( yk(k)*p(k.l) + yk(k+1)*p(k+1,l) )/2; 

end; 
z(l) = H; 
end; 
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% Check the summation of the f(z) density 
F = 0: 

% Change the scale on the z axis scaled In terms of millisecond 
t = zk*1000; 


fork = 1:length(zk)“1. 

F = F + dtz*(z(k) + z(k+1))/2: 

muz(m) = muz(m) + (t(k)+t(k+1 ))*dtz*( z(k) + z(k+1) )/4: 

end; 

F 

muz 

for k = 1 :length(t)-1, 

sigz(m) = sigz(m) + ((t(k)-muz(m))'^2)*dtz*( z(k) + z(k+1) )/2; 
end; 

sigz(m) = sqrt(sigz(m)) 

w = -0.5*( (t-muz(m)*ones(1,length(zk)) )/sigz(m) )^2; 

d_app = (1/(sqrt(2*pi)*sigz(m))rexp(w); 

dtt = t(2H(1); 

tgz(m.:) = z*dtz: 
tgapp(m,:) = d_app*dtt; 

end 

% Output the probability density functions 

figure(1) 

plot(t,tgz/dtt) 

xIabeICCentroid Time of Arrival (milliseconds)’) 
ylabelfProbability Density (1/milliseconds)') 
titleCSampling Frequency 8 MHz') 
text(1.025,20,’3000 Hz') 
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text(1.015,32.'6000 Hz’) 
te)a(1.010.44;9000 Hz’) 
text(1.005.60.’12000 Hz’) 
text(0.96.60.’Peak S/N = 20 dB’) 
axis(I0.95 1.05 0 200.0]) 

% sigzout = sprintf(’%5.4f,sigz(1)): 

% muzout = sprinif(’%5.4f,muz(1)): 

%text(1.25,3,’S/N 15 dB’) 

% text(1.25,2.5,rMean =' muzout]) 

% text(1.25,2.0,[’Std = ’ sigzout]) 

% text(-0.25,5.0,'Calculated Density’) 

% % text(-0.25,3.0,’Gaussian Approximation') 
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% Puldist.m 
% 

% Name: Richard W. Williamson 
% Date; 18 February 1994 
% Date Revised; 18 February 1994 
% 

% This program generates a simulated pulse. The length of the pulse is 
% 1.0 microseconds the maximum amplitude is 1.0, and the sampling rate is 
% 8.0 MHz. The probability density function is calculated for peak signal to 
% noise ratios of 15, 20, 25 and 30 dB. 

% 

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 

% The program outputs one figure. 

% 

% figure(1) A plot of the probability density function for the pulse centroid. 
% 

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 

clear; 

% Initialize the variables 
% 


muz = [000 0]; 
sigz = [000 0]; 

Vk = [0.03162 0.0100 0.003162 0.00100] 

% initialize the variables. 

Dyk = [3;(9-3)/60;9; 3:(9-3)/60;9: 4;(8-4)/60:8; 5;(7-5)/60:7]; 

zk = 3;(7-3)/200:7; 
z = zeros(length(Vk),length(zk)), 

for m = 1 ;length(Vk), 

yk = Dyk(m,;); 
dtz = zk(2)-zk(1); 
dty = yk(2)-yk(1); 


dt = 1.00; 
% 


140 






% The variance of the noise power added to the pulse 
V = \/k(m); 

% The pulse without noise, 
s = [0 0.5 1 1 1 1 1 0.5 01; 

N = length(s); % Length of the pulse. 

% Generation of the noise 
ns = s; 

% Numerator = x Denominator = y; 

% Initialization of the mean and sigma variables 

mux = 0; 

muy = 0; 

sigx = 0; 

sigy = 0; 

cov = 0; 

% Calculate the statistics of the numerator and 
% the denominatior of the centroid equation. 

fork = 1:N, 

mux = mux + k*dt*ns(k); 
sigx s= sigx + ((k*dt)''2)*V: 
muy = muy + ns(k); 
sigy = sigy + V; 
cov = cov + k*dt*\/; 


end; 

sigx = sqrt(sigx); 
sigy = sqrt(sigy): 

rho = cov/(sigx*sigy); % The correlation coefficient 
c_bar = mux/muy; 
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% Calculate the probability distribution for f(zy,y) 
for I = 1:length(zk), 

for k= 1;length(yk), 

y = yk(k); 
z = zk(l): 

g1 = (y*z-mux)/sigx: 
h1 = (y-muy)/sigy; 

fxy = (1/(2*pi*sigx*sigy*sqrt(lH+io'^2)))*... 

exp( (-1/(2*(1-rho''2)) )*(g1''2 -2*rho*g1*h1 + hr2)): 
p(k.l) = fxy; 

end; 

end; 

% Integrate the probability distribution across y to 
% calculate the marginal density for z the centroid variable. 

for I = 1 :length(zk), 

H=:0; 

for k = 1:length(yk)-1, 

H = H + dty*(yk(k)+yk(k+1))*( p(k.l) + p(k+1.l))/4; 
end; 

2 (m,l) = H; 
end; 

% Check the summation of the f(z) density 
% and compute the mean of the distribution 
F = 0; 

% Change the scale on the z axis 
t = (zk-ones(1,length(zk)))/(N-1); 

for k = 1:length(zk)-1, 

F = F + dtz*( z(m,k) + z(m,k+1) )/2; 

muz(m) = muz(m) + (t(k)+t(k+1))''dtz*( z{m,k) + z(m,k+1) )/4; 

end; 
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for k = 1;length(t)-1, 

sigz(m) = sigz(m) + ((t(k)-muz(m))'^ 2 rclt 2 *( 2 (m.k) + 2 (m,k+ 1 ) )/2: 
end; 

sigz(m) = sqrt(sig 2 (m)) 
tg(m.:) = 2 (m,:): 


end; 

figure(1) 

t = (2k-ones(1,length(2k)))/(N-1); 
dtt = t(2)-t(1); 
plot(t,dt 2 *tg/dtt) 

title('Sampling Frequency 8 MHz') 
xlabel('Centroid Time of Arrival (microseconds)') 
ylabel('Probability Density (1/microseconds)') 

% The statistics for 15 dB 
sigzout = sprintf('%5.5f,sigz(1)); 
muzout = sprintf('%5.4f,mu2(1)); 

text(0.545,25.’S/N = 15 dB') 
text(0.545,20,['Mean =' muzoutj) 
text(0.545,15.['Std =' sigzout]) 

% The statistics for 20 dB 
sigzout = sprintf('%5.5f ,sig2(2)); 
muzout = sprintf('%5.4f ,mu2(2)); 

text(0.545,45,'S/N = 20 dB') 
text(0.545,40,['Mean =' muzout]) 
text(0.545,35,['Std =' sigzout]) 

% The statistics for 25 dB 
sigzout = sprintf('%5.5f ,sigz(3)); 
muzout = sprintf('%5.4f ,muz(3)); 


text(0.545.65,'S/N = 25 dB') 
text(0.545,60, [Mean = ' muzout]) 
text(0.545,55,['Std =' sigzout]) 
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% The statistics for 30 dB 
sigzout = sprintf('%5.5f,sig2(4)); 
muzout = sprintf('%5.4f ,muz(4)); 

text(0.545.85,'S/N = 30 dB') 
text(0.545,80,['Mean =' muzout]) 
text(0.545,75,['Std =' sigzout]) 

% axis([0.4 0.6 0 90]) 
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